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SUMMARY 


A method for using system identification techniques to improve airframe finite element models 
using test data has been developed and demonstrated. The method uses linear sensitivity matrices 
to relate changes in selected physical parameters to changes in the total system matrices. 1 he 
values for these physical parameters were determined using constrained optimization with singular 
value decomposition. The method was confirmed using both simple and complex finite element 
models for which pseudo-experimental data was synthesized directly from the fimte element 
model. The method was then applied to a real airframe model which incorporated all ot the 
complexities and details of a large finite element model and for which extensive test data was 
available. The method was shown to work, and the differences between the identified model and 


the measured results were considered satisfactory. 
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NOMENCLATURE 


C 

C A 

Ci 

m 

K 

K a 

ki 

M 

M a 

mi 

q(t) 

Q(t) 

u 

Ui 

Ur 

V 
Vi 
Vr 
W 
Wi 
Wr 

y(0 

Y(t) 

Yi 

y 2 

Z.Z, 

a 

at 

Pi 

Y. 

8 

X 

<P 

<p( f ) 

V 

T| 

<D 

d( ) 


coefficient matrix 
damping matrix 

damping matrix for analytical model 
grouped element damping matrix 

lambda matrix 
stiffness matrix 

stiffness matrix for analytical model 
grouped element stiffness matrix 
mass matrix 

mass matrix for analytical model 
grouped element mass matrix 
displacements at the n degrees of freedom 
n independent forces applied at each DOF 

ZjA 

imaginary component of U 
real component of U 

Z^A 

imaginary component of V 
real component of V 
Z2 

imaginary component of W 

real component of W 

redefined displacement vector 

redefined applied force vector 

-(dUR T M A +dV R T C A +dW R T K A ) 

-(dUi T M A +dVi T C A +dWi T K A ) 

special modal matrices defined in this paper 

eigenvalue 

adjustable physical mass quantities 
adjustable physical damping quantities 
adjustable physical stiffness quantities 
real component of the complex eigenvalue 
eigenvalue 
modal matrix 

imaginary component of the complex eigenvalue 

rth modal column 

n x 2n rectangular modal matrix 

modal coordinates matrix 

modal matrix 

differences between experimental data and analytical data 
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INTRODUCTION 


Background . , . , 

The vast bulk of the work reported to date on identification of structural dynamic systems has 
focused on identifying mathematical models that reproduce test results, but little consideration has 
been given to the physical basis for the identified system equations. Typically, the identification 
procedures make systematic adjustments to the system equation, commonly to the stiffness and/or 
mass matrices but also to the damping matrix, so that the identified eigenvalues and eigenvectors 
reproduce as closely as possible the results measured in tests. The result of this process is almost 
inevitably identified mass, stiffness and damping matrices that are fully populated, that is, which 
have nonzero values for almost all elements. Such matrices, while capable of producing plausible 
eigenvalues and eigenvectors, can nonetheless be physically implausible in the sense that the large 
numbers of nonzero elements throughout the system matrices implies direct connectivity among the 
degrees of freedom that does not exist physically. 

Identified mathematical models that are based on physically implausible system matrices may 
be quite acceptable if the objective of the study is to develop a simulation model. However, such 
results for analysis purposes are generally unsatisfactory because it is difficult or impossible to 
relate specific features of the physical system to the analysis results. This problem is particularly 
troublesome when the objective of the identification of a system model from experimental 
measurements is an accurate system model that, in turn, will be used to make modifications to or 
improvements in the original physical system. Such an example might be the modification of an 
existing aircraft structure to accommodate a new mission. In this case it would be desirable to 
formulate a structural model for the present structure, verify its accuracy against experimental 
measurements, and then use it as the basis for the modifications. When the verification process 
yields identified system matrices that are mathematically acceptable but physically implausible, the 
resulting model may be useless as the basis for future structural modifications. 

The objective of the present work was to develop a method for identifying physically plausible 
finite element system models of airframe structures from test data. The assumed models were 
based on linear elastic behavior with general (nonproportional) damping. Physical plausibility of 
the identified system matrices was insured by restricting the identification process to designated 
physical parameters only and not simply to the elements of the system matrices themselves. For 
example, in a large finite element model the identified parameters might be restricted to the moduli 
for each of the different materials used in the structure. In the case of damping, a restricted set of 
damping values might be assigned to finite elements based on the material type and on the 
fabrication processes used. In this case, different damping values might be associated with 
riveted, bolted and bonded elements. 

The method itself is developed first, and several approaches are outlined for computing the 
identified parameter values. The method is applied first to a simple structure for which the 
"measured" response is actually synthesized from an assumed model. Both stiffness and damping 
parameter values are accurately identified. The true test, however, is the application to a full-scale 
airframe structure. In this case, a NASTRAN model and actual measured modal parameters 
formed the basis for the identification of a restricted set of physically plausible stiffness and 
damping parameters. 

Review of Previous Pertinent Work 

Airframes are generally modelled using powerful finite element analysis packages such as 
NASTRAN that are capable of representing quite detailed aspects of the structural system. The 
accuracy of such models is determined by comparing the analytical results with flight or ground 
vibration test results. In the case of helicopter airframes, several recent efforts have focused on the 
correlation of NASTRAN model data with ground vibration test data 1 ' 3 . The conclusions reached 
in these studies suggest that in cases where there is some degree of correlation, the model 
frequencies compare favorably with test frequencies, but generally only in the low frequency range 
below about 15 Hz u2 . The frequency response functions at selected locations also compare 
reasonably well in this range. Outside this range the comparisons are generally unsatisfactory, and 
the eigenvectors do not usually compare favorably in either range. 

Although there have been numerous contributions to the literature in the area of the 
identification of structural dynamic systems 4 * 25 , the majority of reported methods are based on 
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simply adjusting the elements of one or more of the K, M, and C matrices. While this approach is 
capable of yielding a system matrix whose eigenvalues and eigenvectors suitably match measured 
results, the methods generally lose all physical interpretability inherent in the original K, M and C 
matrices by not maintaining relationships among elements dictated by the model topology. These 
difficulties are compounded for large-scale models with thousands of degrees of freedom. 

The reported papers, that address this problem are by Meirovitch and Norris 25 ' 27 , Lim 28-30 , 
Hajela 31 , Zimmerman 32 , Hickman 33 et.al. and Chen and Garba 33 . In reference 25, the problem 
of parameter identification in distributed parameter system has been addressed. It is assumed that 
the response of a structure is measured and the physical parameters of the system are identified by 
using least square principle. In reference 26, a perturbation technique has been discussed for 
parameter identification. As a basis for the perturbation procedure, it has been assumed that there is 
a prior knowledge of approximate values of system parameters. The developed perturbation 
technique depends on measured system response to a harmonic excitation in a frequency domain. 
Again, the identification is based on least square principle. In reference 27, the method has been 
extended to a Raleigh-Ritz type of method where the parameters have been assumed to be known 
functions with undetermined coefficients. The developed theory considers linear viscous damping. 
However, the numerical examples do not include damping and are based on simulated 
measurements that have been obtained by numerical methods. Lim has used a submatrix approach 
to update stiffness matrices in reference 27. In reference 28, both stiffness and mass matrices have 
been updated simultaneously. In reference 27 and 28, Lim has attempted to retain the physical 
significance of the stiffness and mass matrices by use of sub-matrices that represent the stiffness or 
mass matrices of individual or group of elements. In reference 29, Lim has used a method that is 
similar to that of references 27 and 28 to identify damages to structures. In all these works, Lim 
has considered undamped system matrices and has used numerically simulated eigen values and 
eigen vectors to validate the theory. In references 31-33, the identified physical parameters are the 
structural failures. These failures have been related to changes in stiffness matrices. The algorithm 
is based on partial inverses and optimum least squares expansion. Damages are identified only by 
observation of global changes in properties. Quantities such as parameter sensitivity or parameter 
perturbations or methods such as singular value decomposition or constraints to limit parameters to 
physical space have not been used. Also, these methods have not been applied to large scale 
system matrices. In reference 32, however, a pattern recognition technique has been used on the 
basis of measured response of the structure. 

In our work, we would like to consider measured frequencies and modes as inputs to identify 
system parameters that can be related to physical variables. Our objective is to identify large scale 
systems like helicopters that have been specified by a known finite element model. Because of the 
sensitivity of errors in experimentally determined modes to system parameters, it is also necessary 
to impose constraints on identified parameters to occupy specified parameter space. For example, 
we can not have a negative modulus or a negative damping in a passive structural dynamic system. 
In using measured data, we found that the identification process needed the use of singular value 
decomposition methods. Because we are using experimentally generated data, the examples 
include general linear damping matrices. No restriction of proportional damping has been 
imposed. 

Kuo and Wada 35 used nonlinear sensitivity coefficients (NSC) in the identification procedure. 
Their sensitivity coefficients are between the system parameters and eigenvalues. In the present 
work the interest is in the change of system matrices as a function of physical variables of the 
structure. A different type of sensitivity coefficient between system matrices and physical variables 
has therefore been developed. 

The most significant achievement in the present work 40 is to preserve the physical 
interpretability of the M, C, K matrices so that the identification can provide evidence of possible 
sources of erroneous modeling and point to specific regions of the model that are unduly sensitive 
and need additional consideration in modeling. The identification procedure developed in this paper 
is capable of adjusting physical quantities such as boundary conditions, moments of inertia, 
stiffnesses, damping or other selected physical parameters. 

Mathematical Model 
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Basic Equations 

Any linearly elastic structural system with n discrete degrees of freedom and with general 
viscous damping (either proportional or nonproportional) can be represented by n coupled ordinary 
differential equations that can be written in the following form 37 : 

Mq(t)+Cq(t)+Kq(t)=Q(t) (D 

where M, C, and K are symmetric n x n inertia, damping, and stiffness matrices, respectively. In 
this formulation, q(t) are the displacements at the n degrees of freedom and Q(t) are the n 
independent forces applied at each degree of freedom. 

In the case of undamped or proportionally damped systems, there are n complex conjugate 
pairs of eigenvalues and n distinct modes which are orthogonal with respect to M and K. Using a 
transformation matrix of the form: 


qCO^TlCt) 

will allow decomposition of the original system equations (Eq. 1) into n decoupled equations that 
are straightforward to solve. 

This transformation cannot be applied to the general nonproportionally damped problem in the 
same manner because for this case there are 2n complex modes, <p^, and consequently 2n modal 


coordinates, t| r (t), but there are only n physical coordinates, qi(t). 

One can overcome this difficulty by writing Eq. (1) as a set of 2n ordinary differential 
equations in the form: 



( 2 ) 


If one then defines: y(t)= { } and Y(t)= { } , the above equations can be written as a set 

of 2n first order ordinary differential equations: 


rO Mi 

-M 

0 l 

Lmc>H 

. 0 

K J’ 


(3) 


This formulation has the advantage that the modes obtained from the solution of the homogeneous 
equations, obtained by letting Y(t)=0 in Eqs. (3), are orthogonal, and hence can be used in 
conjunction with the expansion theorem to obtain the solution of the nonhomogeneous problem. 
The solution of the homogeneous equations is obtained by assuming as before a solution in the 
form: 


y(t)=3>e at 


(4) 


where O represents the spatial component of the solution and is a vector consisting of 2n constant 
elements. The corresponding eigenvalue problem can be written as: 


0 Ml r-M 0 i f01 

“LmcH# kJ'Ho) 


(5) 


The solution of the eigenvalue problem yields 2n eigenvalues, otr, and 2n eigenvectors 


0>T = 




r=l,2,...2n. 


( 6 ) 


Equations (5) and (6) provide the solution to Eq. (1), but in order to simplify the computational 
work, it is convenient to formally separate these complex equations into real and imaginary pairs. 
Following the approach introduced by Cheng 10 , the real and imaginary components of the 
eigenvalues and eigenvectors are defined, respectively, as: 


^r^Rr+j^If 
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and in addition the following modal matrices are defined: 

Zi = (Re(oti4>i), Im(ai<t>i), Re(a2<t>2), Im(a2<t>2). •• 

Re(a n <M. Im(a n 4>n)) 

Z2 = (Re(4>i), Im(<l)i), Re(<|>2), Im(<J>2), ... 

Re(<t> n ). Im(<l>n)) 

Then the new system equation, Eq. (5), can now be rewritten with purely real terms in the form: 


( 8 ) 


(9) 


( 10 ) 


4mcH£Mo M kH%H»} 

where the eigenvector matrix. A, is a block diagonal matrix with blocks 
Sr 

.-a Sr _ 

along the diagonal and zeros elsewhere. Equation (9) can be further simplified by the introduction 
of U, V, and W as follows: 

MU+CV+KW=0 (11) 

where U = ZiA, V = Z2A and W = Z2 or expUcitly: 

U = (ReCai 2 ^!), Im(oti 2 <J>i), Reta^te), Im(ct 2 2 <l> 2 )* 

... Re(a n 2 <t> n ), ImCctn 2 ^)) 

V= (Re(ai <}>i ), Im(ai<()i), Re(a 2 <fe), Im(a 2 <J> 2 ), ... 

Re( 0 Cn <t>n), Im(a n <l)n)) 

W = (Re(<J»i), Im(<t»i ), Re(<f> 2 ), Im(<}> 2 ), ... Re(<I> n ), 

Im«M) ^ 12) 

Finally, Eq. (1 1) can be separated into explicit real and imaginary equations in the form of the 
following two equations. 

MU r +CV r +KWr=0 (13) 

MUi+CV!+KWi=0 (14) 

These equations are same as Eqs. (5), but they do not include complex variables. For the 
identification procedure, it is much easier to use these equations than to use Eqs. (5) directly. 

lde To beginTuppose ^at the mass, damping and stiffness matrices for the initial analytical model 
are given by M A C A and K A , respectively, and the identified mass, damping and stiffness matrices 
are given by M, C and K. In a similar manner, the eigenvectors and eigenvalues for the analytical 
model are given by U A , V A and W A , while U E , Ve and W E are the eigenvectors and eigenvalues 
determined from test data. From these definitions it follows that the relationship between the 
identified model (based on the test data) and the analytical model can be written as: 


M=M A +dM, C=C A +dC, K=K A +dK 


(15) 


U E =U A +dU, V E =V A +dV, W E =W A +dW ( 16 ) 

where dM, dC, dK, dU, dV and dW are the changes. The identified model satisfies equations 
(13) and (14), so substituting equations (15) and (16) into equation (13) and (14), yields. 
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( 17 ) 


dU R T M A +dV R T C A +dW R T K A = 

-(UER T ,VE R T ,W ER T )(dM,dC,dK) T 


dUi T M A +dVi T C A +dWi T K A = 
*(Uei t »V ei t »W Ei T )(dM,dC,dK) T 


These equations can be combined into 
[" Uer 7 Ver 7 Wer t | 

Luei 7 Ve! 7 Wei 7 J1 


the following form: 



(18) 

(19) 


where 


Y! =-(dU R T M A +d V r t C A +dW r t K a ) 

Y 2 =-(dUi T M A +d V i T C A +d W i t K A ). ' J 

The right side of these equations is known, since Ma, Ca* and are given by the analytical 
model and dU R 7 , dV R 7 , dW R 7 , dUj 7 , dVj 7 , and dW] 7 , which are the differences of the 
eigenvalues and eigenvectors between the analytical model and the experimental data, are known. 
Finally, the matrix 


'U ER 7 Ver 7 Wer 7 ' 
.Uei 7 V m T Wei 7 . 


contains only experiment data. 

The solution to these equations are the changes of dM, dC and dK. Because of matrix 
symmetry, the number of unknowns in Eq. (19) is 3 n(n+l)/2. The number of equations depends 
on the number of known experimental modes. Suppose this number is m, then the number of 
equations are m x n. If the number of the equations is larger than or equal to the number of 
unknowns and the rank of this matrix is equal to 3 n(n +l)/2, normal least square methods can be 
used to solve these equations. Otherwise, singular value decomposition, or constrained 
optimization can be used to solve Eq. (19) for the changes dM, dC and dK, and these results can 
then be substituted into Eq. (15) to determine the identified M, C and K matrices. It should be 
noted that this approach is capable of handling nonproportional damping and underdetermined 
problems in which fewer modes are measured than are computed from the analytical model. 

At this stage the usual identification procedure can be performed. The values of M, C and K 
can be put into the system equation, Eq. (1), and the experimental data can then be reproduced. 
However the identified M, C and K cannot be related to particular physical quantities in the actual 
airframe, because the changes occur throughout the entire M, C and K matrices . In order to 
preserve the physical interpretability of the identified system, it is necessary to develop a 
relationship between dM, dC and dK and adjustable physical quantities such as boundary 
conditions, moments of inertia, stiffnesses or other selected physical parameters. To this end, 
assume that each of the system matrices can be decomposed into the form: 

N m N c N k 

M^Tm.a,, C=X c i(5i, and K=X k iY> ( 21 ) 

i=l i=l i=l 

where ctj, pi and Yi are adjustable physical quantities and mj, q and kj are grouped element 
matrices with common physical quantities. 

For example, in the finite element model of actual airframe, there is an ej-th element , (see Fig. 
1). The portion of the stiffness matrix that describes bending in the xz plane of an element, 
assumed to be a principal plane (Fig. 2), in NASTRAN, is given by 
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R 


-^R -R -jR 
l2 r; + EI V Ir — r ^ 

4 r+ j 2 r 2 r- j 

sym R ^R 

J 



where R - (j^G+ i2EI y ) * 
physical quantity, yk, then 


If the modulus of elasticity, E, is taken here as an adjustable 
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K ej = 


R 


L sym 
r 


-^R -R -^R “* 

i 2 r Ely I R 1 ? r eIy 

4 r + j 2 R 2 R 1 

R iR 

1 2 „ Ely 

1 n 


1 

‘2* 


-r 


l 2 Iy 1 l 2 

4 t+ 1 f I*' 1 


L s y ra 

f2(l+v)l l3 


1 

• f 

l 2 Iy 
4^+f -J 


7k = k«j7k 


(23) 


where Yk= E and r = . Suppose there are n elements which have the same E 

so that it is possible to express the stiffness as: 
kk = Ikej 

Cj=l 

When the modulus changes from E to E+dE, the corresponding change in 7k is to 7k+ d 7k. 
Considering all different 7k, K changes from K to K+dK where 
N k 

dK=£Mu 

k=l 

Similar procedures can be generalized to include the damping, other stiffness parameters, and 
mass. 


N m N m 

dM=V ^da^X^ 0 ^ 

jL^dOi i=i 


i=l 

N. 


N„ 


dC=\^ |^dPi=X c i d Pi 
£^ 0 $ i i=l 


(24) 


i=l 

N, 


N,. 


dK=V ~d7i = Xki d 7i 

i=i 

i=i 


Substituting these into Eq. (19) yields a set of linear algebra equations with unknowns doq, dpi 
and dy»: 
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UrT 


..VrT^-...WrT^. 

api ayi 


9M 
3ai 

'8m ,,_ T 9c w 

api ' an 


UlT^...ViT- 

9a i 


r dai 


daN m 

dpi 


>= 


dPN c 

dyi 


(25) 


v. dyNk ^ 

The number of unknowns in this equation is much less than the number of unknowns in Eq. (19), 
and also all the unknowns in this equation have physical meaning in the real structure. 

However, neither Eq. (19) or Eq. (25) can be solved directly since the numbers of unknowns 
and equations are not equal in most of the cases. There exists a number of techniques for dealing 
with sets of equations that are under or over-determined or with matrices that are either singular or 
else poorly conditioned. The singular value decomposition, or SVD method 36 , is one of the most 
powerful ways to handle these problems. In the present study it was employed to compute 
solutions to Eq's. (19) and (25) which are highly under-determined for most practical situations. 

In this case the SVD method provides a least square type of solution to the problem. 

In most cases, the selected physical parameters must also be restricted to positive values in 
order to make sense physically. However, the identification procedure outlined above cannot 
guarantee that the identified values will all be positive. This is of particular concern when the 
parameters are proportional to mass, an elastic modulus or a damping coefficient, all of which must 
be positive for the systems typically considered. Using a constrained optimization method, this 
problem can be eliminated. The present problem can be posed as one of minimizing 

f=dai+da2+-+daN m +dpl+..+dPN c + dYl + ” + dYNk ( 26) 

subject to the constraints 

r dai "I 


UrT 


Ui 


^...V R T^...W R TfS. 
9ai 9Pl 9yi 

tM Vi t — WtT^K 


9ai 


— ... Wi J 

9Pl 9yi 


daNm 

dpi 


dPN c 

dyi 





Y2 


v dyNk J 

and 
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d«i>0, da 2 >0, ..da Nm >0, dpi>0,..dp Nc >0, dyi>0, ..dyN k >0 

The feasible solution (dc^, da 2 , .... da Nm , d$ u .... dp Nc , dyi, .... dy Nk ) to this problem yields 
the identified selected physical parameters. 
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APPLICATIONS 


System Identificatio n Procedures . , , „ . ., 

The method described above was applied to several practical examples. For these cases, the 
analytical finite element model for the structures was assumed correct and was developed using the 
NASTRAN finite element analysis package 38 . Then, values for selected physical parameters in the 
model were identified on the basis of measured experimental data (eigenvalues and eigenvectors) 
so that the analytical model more accurately represented the real structure. The assumption for this 
procedure was that the identification process could be applied in an iterative fashion by making 
successive small modifications to the selected physical parameters until satisfactory agreement with 
experimental results was obtained. For the i-th iteration, there are the following relationships. 

Mi = M^ 1 + dM, C» = C*' 1 + dC, Ki = K^ 1 +dK 

ui = U 1 ' 1 + dU, yi = yi- 1 + dV, wi = W*' 1 +dW, (27) 

Substitute these into equation (25), we can obtain 

f dal 1 ^ 


Ur 


1»! vr^...wr^ 


9ai 


api 


ayi 


... VlT^l ... 


9oq 


dpi 


dyi 



(28) 


V dYNk 1 ' 

where 

Yi' = -(dU^Mi-i + dv£ci-l + dW^Ki* 1 ) 

Y2 -(dU^Mi- 1 +dv’[c i - 1 + dW^' 1 ) (29) 

The convergence criteria was formulated as follows: 

(1) Check the physical parameter differences doci, dcx2, ••» dccNm, dPi, .., d(3Nc» dyi, dyNk 
either manually or programmatically. If these physical parameter differences are smaller than a 
tolerance value, the identified physical parameters are obtained. 

(2) Check the differences of the experimental eigenvalues and the i-th iteration analytical results 
which are obtained after running NASTRAN. If the differences are smaller than a tolerance 
value, the identified system is obtained. 


Simple Numerical E xample . , , . . , 

In order to verify the proposed approach, the identification procedure developed above was 
applied first to a very simple finite element model with only a few degrees of freedom. It is a 
simple variable cross section straight rod with fixed ends, and it contains all the desired parameters 
to be identified such as mass, stiffness and damping. It was modeled using 9 rod elements with 
lumped masses at each node as shown in Fig. 3, and representative values were assumed for all 
elements and mass properties. For the purpose of defining the damping, the elements were 
segregated into 4 groups and a different damping coefficient was specified for each group. 
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Figure 3. Simple Numerical Example 

The assumed physical properties were defined to be typical of an aluminum rod. The length 
was 228.6 mm (9") and the modulus E =68.9 GPa (10e6 psi). The concentrated masses at each 
node were those given in Table 1. This model was then employed to generate eigenvalues and 
eigenvectors which were used to represent “measured” test results. In all the cases presented in 
this report, the calculations were done in US customary units and the results converted to SI units. 
As a result, some of the percentage figures may be slightly in error due to numerical roundoff in 
the conversion of units. 


TABLE 1. 

ASSUMED PHYSICAL PROPERTIES FOR SIMPLE NUMERICAL EXAMPLE 


Index 

Mass at 
Node 
(kg) 

Element 

Stiffness 

(MN/m) 

Element 
Damping 
Coefficient 
(kN s/m) 

1 

21.89 

367.7 1 

28.34 

2 

22.76 

369.5 ] 


3 

23.64 

lau 

28.34 S 

4 

24.51 

373.0 


5 

25.39 

£2Z2MI 

27.58 

6 

26.26 

376.5 

27.14 

7 

27.14 

378.3 1 

1714 “ 

8 

28.01 

HH 

^ H 

9 

28.89 

381.8 

26.27 

10 

29.76 




The objective of the identification was to determine the physical parameters such as damping 
constants (c or Q and the cross section area for each rod element There are three different cases to 
start to consider with this system. In the first case, the mass and the stiffness matrices were 
assumed to be accurate, and four damping parameters were identified assuming zero as initial 
analytical values of the damping matrix. The identified damping parameters are listed in Table 2 . 
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TABLE 2. 

CASE I: IDENTIFYING THE DAMPING PARAMETERS 


Damping 

Parameter 

Exact 

Value 

Initial 

Value 

Identified 

Value 

Error 

(%) 

ci 

28.34 ] 

HK9I 


msm 

C 2 


WMSMM 

wi-xma 

H 

C3 

mmm 

IH&H 


W 

C4 


WMSSM 

fascia 

mm 

C5 




mat i 

C 6 

27.14 

HKSlH 

27.1431 

w 

ci 

27.14 

bb 


SB 

C 8 

26.27 

la 


-0.0005 

C9 

26.27 


\ 26.2676 

-0.0005 


The second case was to identify the stiffness parameters assuming accurate values of mass and 
damping parameters which were the same for all elements ki = 376.5 MN/m. The identified 
stiffness parameters are listed below. 


TABLE 3. 

CASE II: IDENTIFYING STIFFNESS PARAMETERS 


Stiffness 

Parameter 

Exact 

(MN/m) 

Initial 

(MN/m) 

Identified 

(MN/m) 

Error 

(%) 

kl 

367.7 

376.5 1 

367.7 

0.0000 

k2 

369.5 

376.5 

369.5 

-0.0005 

k3 

371.$ 

376.5 

371.3 j 

0.0000 

k4 

373.0 

376.5 

373.0 

0.0000 

k5 

374.8 

376.5 

374.& 

0.0000 

k6 

” 376.5 1 

376.5“ 

376.5 1 

0.0000 

ki 

378.3 

376.5 

378.3 

0.0000 

k8 

380.0 

376.5 

380.0 

0.0060 

k9 

$si.r“ 

376.5 

381.8 

0.0000 


In the third case, both the damping and the stiffness parameters were identified under the 
assumption of accurate mass value alone. The elements of the initial damping matrix were assumed 
to be zero, and the stiffness parameters were assumed to be the same for all elements (ki- 51o.d 
MN/m). The identified damping and stiffness parameters are listed in Tables 4 and 5. 
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TABLE 4. 

CASE m(A): DAMPING PARAMETERS 


Damping 

Parameter 

Exact 

Initial 

| 

Error 

(%) 

ci 

28.34 



mm 

C2 

28.34 

0. 


0.0000 

C3 



28.34 

0.0000 

C4 

27.58 




C5 

27.58 | 

o: 

i mm 

[iMGI 

C6 

esbsm 

HH9H 

EH ma 

rnomia 

C7 

mmm 

HK9H 


rnriTiiEi 

C8 

26.27 

HBSSHI 


BOBBISH 

C9 

WSS5M 

^H£HI 

PHEK 

OTtHIttl 


TABLE 5. 

CASE m(B): STIFFNESS PARAMETERS 


Stiffness 

Parameter 

Exact 

(MN/m 


Identified 

(MN/m) 

Error 

(%) 

kl 

367.7 

376.5 

mr&Mcm 

KSU 

k2 


bees 


msm 

k3 

■ana 

Eras 


mm 

k4 

waw 

ms&a 

373.044 

mm 

k5 

374.8 

376.5 

muLMita 

msm 

k6 

376.5 


376.521 

mm 

kj 

378.3 

mu.*m 

■twawflil 

n 

k8 

380.0 



BlIlIUI 

k 9 


376.5 


mm 


All the results were obtained after only one iteration. For these simple cases the method 
accurately identified the selected physical parameter values (damping and cross section areas). 

A pplication to AH-1G Model 

A NASTRAN finite element model (FEM) for the AH-1G helicopter airframe has existed for a 
long time and was originally developed by Bell Helicopter Textron Inc. It is basically composed of 
two parts, one is stiffness modeling for idealizing the structures and the other is weight modeling 
for distributing weights to grid points. There are 4405 different elements with a total of 2764 
degrees of freedom in the basic full model. A reduced model, based on Guyan reduction, contains 
only a total of 63 physical degrees of freedom. 

Normally, the input and output data files from NASTRAN dynamic analyses are specially 
formatted and are quite large for a large finite element model such as the full AH-1G model. For 
convenience and accuracy, the present system identification programs were designed to 
automatically read NASTRAN output files and create NASTRAN input data deck files. At each 
step in the iterative identification procedure, the new modified physical parameters were put into 
the NASTRAN model bulk data in order to generate the required analytical results, such as 
eigenvalues, eigenvectors and other parameters, for the next iteration. 

The mass, stiffness and damping matrices defined with respect to the internal degrees of 
freedom are not normal NASTRAN output data. However, such results can be developed by 
using appropriate Direct Matrix Abstraction Programming (DMAP) utilities so that the identification 
program can automatically get this NASTRAN output data (see Appendix B). 
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Results Using Simulated Test Data . Wq1 - 

The NASTRAN model of an AH-1G airframe includes 4405 different elements with a total of 
2764 degrees of freedom. In order to make sure that the identification procedure was appropriate 
to a such big model, the use simulation has been chosen to begin with. For this identificauon, the 
mass and stiffness properties of the analytical model were considered to be accurate, and 
nonproportional damping properties were identified. The physical damping parameters were 
associated with 8 distinctly different types of materials and structural fabrication techniques used in 
the airframe (e.g., aluminum, steel, riveted, welded, bolted, etc.) and one of these damping values 

was associated with each of the model elements using Eq. (1 1). 

For this case, the test data were synthesized from the original NASTRAN model assuming 
small values for the extension and rotation viscous damping coefficients (kN-s/m and N-s/rad 
units): 

TABLE 6 . ASSUMED INITIAL PHYSICAL DAMPING VALUES 


Extension 

Rotation 

Cl =5.2531 

C 5 =93.4 

C 2 = 8.756 ! 

C 6 = 155.7 

C 3 = 1.751 

C 7 = 31.14 

C 4 = 1.226 

C 8 =21.80 


The synthesized data included 24 modes of which 6 were rigid body modes, and the frequency 
range was from 0.0 to 30.2 Hz. The dimension of the mass, stiffness and damping matrices was 
2764 x 2764. The initial values of the physical damping parameters for the analytical NAS 1 KAN 
model were taken to be zero, and the results for the identified values are shown below: 

TABLE 7. IDENTIFIED PHYSICAL DAMPING PARAMETERS 


Parameter 

Initial 

Identified 

£1 

5.253 j 

5.429 

C'2 

8.756 

RT490 1 

£3 

1.751 

r “ 1.746 

C4 

1.226 1 

~ 6.069 

C5 

“ 93.4 

55.91 

C6 

155.7 

13035 

C 7 

31.14 

653)6 

C8 

21.80 

55.91 


The error in the identified damping parameters as a function of the number of matrix elements for 
each of the 8 damping types is shown in Fig. 4. The error for those element types with more than 
100 elements present is quite low, but it is much larger for those types with only a few elements 
present in the complete finite element model. The largest error was associated with what appeared 
to be elastomeric materials. 
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Fig. 4. Error in Damping Estimate as a Function of Number of Matrix Elements in Model 

This simulation confirms the identification procedure for a complicated but yet well-defined 
example. If the assumptions such as nonproportional damping are correct for the airframe and if 
the experimental data are of high quality, the physical damping parameters can be identified from 
the test data. 

Actual AH-1G Data 

Actual test data for an AH-1G airframe were provided by Bell Helicopter Textron Inc. based on 
ground vibration tests and included both resonance dwell and FRF (frequency response function) 
data. The experimental data were available for 8 different configurations of the AH- 1G that were 
tested. The principal difference between the tests concerned the degree of complexity of the actual 
airframe tested. At one extreme, the bare airframe without most attachments was tested while at the 
other extreme the complete airframe with all attached mechanical components was tested. The test 
data from the most complex airframe configuration (with all difficult components present) showed 
the poorest agreement with the corresponding analytical model, while the data from the simplest 
test airframe showed the best agreement. 

For this study, the test data from the most complex airframe configuration was used. Only the 
FRF data were employed, and the complex eigenvalues and eigenvectors for some 7 triaxial modes 
were obtained from the FRF data by using the TDAS® curvefitting program 39 . The experimental 
data were provided as FRF's in TDAS universal file format, and the results generated by TDAS 
were complex eigenvalues and eigenvectors. 

Before use in the identification program, the eigenvectors were normalized. Two options were 
used to normalize both experimental and analytical eigenvectors. One was to normalize the 
eigenvectors to the same point, and the other was to normalize based on the minimum deviation 
between the analytical and experimental eigenvectors. 

The full finite element model for AH-1G airframe, as mentioned in the previous section, has a 
total of 2764 degrees of freedom which is very large for the identification procedure. In order to 
keep the problem tractable, a Guyan reduction was used in the present application to reduce the 
analytical model to a total of 63 degrees of freedom, which corresponded to the 23 distinct 
locations on the airframe at which experimental measurement were made. The error due to the 
reduction in degrees of freedom from 2764 to 63 is shown in Table 8. 

TABLE 8. 

EIGENVALUES (FREQUENCY) (WITHOUT ANY DAMPING) 


Test 

Full 

Model 

Error (%) 

Reduced 

Model 

Error 

(%) 

7.2475 

7.6734 

5.877 

7.6932 

6.150 

8.0458 

8.3467 

3.740 

8.4026 

4.435 

15.9539 

14.6722 

-8.034 

15.825 

- 0.6 10 

17.2174 

17.3701 

0.887 

17.7&4 

3.294 

23.7396 

20.7955 

-12.392 

22.881 

-3.606 

24.6675 

25.7955 

4.573 

28.238 

14.475 

32.6848 

31.7526 

-2.852 

33.786 

3.369 


Initially, both the full and the reduced models were used as analytical models. Using the actual 
experimental data, the physical parameters in the analytical models were obtained using the 
iterative procedure outlined earlier. The initial results for both the full model and the reduced 
model included several negative identified damping parameters which were obtained using the 
singular value decomposition method when either zero or positive initial guess values were 
assumed for the analytical model. Physically of course, the damping parameters should be greater 


® TDAS (Test Data Analysis) is a part of I-DEAS which is a computer-aided engineering product 
of Structural Dynamical Research Corporation (SDRC). 
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than zero, but mathematically, the identification procedure is oblivious to this constraint. The 
constrained optimization procedure outlined earlier was therefore used in order to overcome this 
problem. In addition, the reduced model was used in most of the identification cases, expect when 
otherwise stated, because of the small error and big savings in computational time. 

The complete system identification was carried out in two steps. The first step was to identify 
the stiffness, and for this process the initial damping values were assumed to be zero. The second 
step was to use the stiffness values obtained from the first step to identify the damping values. 

This was done under the assumption that the greatest change in natural frequency can be obtained 
by changing the stiffness parameter, while changes in the damping parameters will only fine-tune 
the eigenvalues but will obtain accurate modal damping estimates for the structure. 

At the first step, four stiffness parameters associated with elastic moduli for four principal 
materials used in the airframe were selected to be identified. After two iterations, the differences 
between the identified and the initial moduli and the analytical and experimental eigenvalues were 
those shown in the following tables: 


TABLE 9 

IDENTIFIED MODULUS VALUES 



Initial 

(GPa) 

After 

first 

iteration 

(GPa) 

Change 
from initial 
value (%) 

After 

second 

iteration 

(GPa) 

Change 
from initial 
value(%) 

mat.-l 

22.1 

21.7 

-1.81 


23.03 


72.4 

M&M 

0.19 

9.417 



■iffijTtl 

190.5 

-4.75 

28.435 

IGLSUH 

IhENB;! 

■MtM 

112.2 

-7.01 

1 

i 



TABLE 10. 

IDENTIFIED EIGENVALUES (FREQUENCY) 


Test 

Original 

Error 

(%) 

After 

first 

iteration 

Error 

(%) 

After 

second 

iteration 

Error 

(%) 

7.247 

7.693 

6.15 

7.686 

6.05 

■rIEEfiM 

2.47 

8.046 

8.403 

msm 


MXtm 



15.95 

15.82 

WMM 

15.795 

E!£21 


ME$m 

17.22 

17.78 

wmm 

17.762 

Ha 

17.180 

-0.22 1 

23.74 


-3.61 

mwkyrjm 


21.819 


24.67 

28.24 i 

14.5 


14.3 


11.6 

32.68 

33.79 

3.37 1 

| 33.805 

3.43 

32.481 

-0.62 


As the second step, the damping parameters were identified for the previously identified 
stiffness conditions. Initial estimates for the damping parameters were developed by assuming a 

nominal damping ratio, i^=5%. For the extensional elements, it was therefore assumed that the 
initial viscous damping values would be ce=17.5 for all extensional viscous damping, and that for 
the rotational elements (assuming the cross section area to be a circle) it would be cr=222 for all 
rotational damping. 

After one iteration, the results shown in Table 1 1 were obtained. 
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TABLE 11. 

FINAL RESULTS FOR AH-1G MODEL 


Mode 

Test 

||§|§i 

NASTRAN (final) 



C(%) 

I9SSI 

1S8I 

C(%) 

First Lat Bending 

7.247 

2.19 

wkm m 

7.425 

3.00 

First Vert Bending 


1.56 


8.057 

4.55 

Second Lat Bending 

15.954 


BHiM 


■BS 

Second Vert Bending 


mssm 

mm 

■im 

7.48 i 

Fuselage Torsion 


1.70 

KM! 

KliM 



24.667 

msam 



6.25 f 


32.685 

1.95 

33.74 


mi&JUm 
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CONCLUSIONS AND RECOMMENDATIONS 

A structural dynamic system identification procedure that is capable of identifying physical 
parameter changes has been developed. The changes in physical parameters of the system can 
therefore be related to observed experimental data. In the examples considered, physical 
parameters, such as the damping constant of a material that will result in a nonproportionally 
damped system, the modulus of elasticity of a material, and the dynamic stiffness of a beam 
element have been identified by using the experimentally obtained frequency response functions, 

modes and eigenvalues. . . . . 

Following the validation of the developed procedures by using synthesized data on a small 
model, the method was applied to a large-scale NASTRAN finite element model of a helicopter 
airframe. Both synthesized data and observed experimentally identified modal data were used. 
Again, modulus of elasticity, stiffness and damping constants were the parameters considered for 
the four representative materials used in the airframe. With the exception of one material that a 
been used to construct a very small number of components, other material constants were identified 
reasonably accurately where synthesized data were used. When experimental modal data were 
used, the modal parameters calculated from the identified model did not yield the experimentally 
observed modes only in cases where the initial a priori finite element model output and the 
experimental model output differed considerably. When experimental output and the a pnon model 

output were reasonably close, the results of the identification were satisfactory. 

Even though the method was shown to work and the difference between the identified model 
and the experimental observations were considered satisfactory in some cases, there are some other 
cases that need improvement to make the procedure applicable to a structural dynamic design 
nrocpss' 

( 1) While the numerical processes were improved and refined, no similar improvements in the 
quality of the test data could be realized. One result of this problem was that it was relatively 
difficult to match measured eigenvalues and eigenvectors with corresponding analytical 
values. Quite often, the measured and initial eigenvalues matched closely while the 
eigenvectors differed considerably, and the identified eigenvectors were not significantly 
closer in agreement. For this reason it is necessary to consider other experimental data, such 
as the AH- 1G dwell data, which have been acquired by other means. 

(2) In cases where selected portions of experimental data and a priori analytical data oilier 
significantly while a large amount of experimental and analytical data are close together, it is 
necessary to minimize first the large errors by using H« type of identification before using the 
least square analysis with singular value decomposition. 

(3) it is important that a larger group of identifiable parameters be considered. 

(4) It is necessary that we examine the convergence and accuracy of the complete process. 

(5) We have used linear sensitivity coefficients. Accuracy and convergence may require nonlinear 

sensitivity coefficients. 

(6) The real damping in a structural dynamic system may not be linear viscous damping with a 
nonproportional behavior. It is necessary to include other types of damping mechanisms. 

(7) As pointed out by Bell’s DAMVIBS conclusions 1 , nonlinearity is important m considering 
selected components of the airframe. 

(8) We should also examine the experimental parameter estimation processes used to determine 
modal parameters used as inputs to the identification process. 
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APPENDIX A 

SYSTEM IDENTIFICATION PROGRAM LISTING 

The following pages contain a listing of the current version of the program used to carry out the 
structural system identification described in this report. The program is written in the CDC version 
of Fortran 77 and was run on a CDC Cyber 180-990 running under the NOS/VE operating system. 
The program requires the use of the IMSL Scientific Library (Version 1 1) in order to carry out the 
singular value decomposition and the constrained optimization procedures. 

The program was used with Version 66C of MSC/NASTRAN which was also run on the same 
computer system and was used to solve the structural dynamic eigenvalue problems. 
MSC/NASTRAN was used to run the initial eigenvalue problem and all subsequent iterative 
solutions. As a result, the program listing also includes the necessary I/O calls needed to operate 
directly with MSC/NASTRAN input and output data files. The program also requires the 
experimentally determined eigenvalues and eigenvectors to be present is separate input files for 
each eigenvalue. 

A. Summary of Parameters Used bv th e Program 


Parameter 

Dimension 

LAMBDR, LAMBDI 

MD 

LAMBDTR, LAMB DTI 

MD 

WAR, WAI 

NDxMD 

WTR, WTI 

NDxMD 

VTR.VTI 

NDxMD 

DUR, DUI 

NDxMD 

DVR, DVI 

NDxMD 

DWR.DWI 

NDxMD 

MASS 

NV 

STIFF 

NV 

DAMP 

NV 

DC, D^ 

NV 

DK 

NV 

COEFF 

NMxID 

Y 

NM 

WK 

ID2 

BETA 

ID 

GAMMA 

ID 

MODES 

22x2 

NTESTk (k= 1,2,5) 

20x4 

NT 

7 

IRTYPE 

12 


BL,BU 12 

A 12x12 


Definition 

real and imaginary parts of analytical eigenvalues ^a 

real and imaginary parts of test eigenvalues X 

real and imaginary parts of analytical eigenvectors Wa 

real and imaginary parts of test eigenvectors W 

real and imaginary parts of Va 

real and imaginary parts of U-Ua 

real and imaginary parts of V-Va 

real and imaginary parts of W-W a 

mass matrix 

stiffness matrix 

damping matrix 

Ci matrix 

kj matrix 

coefficient matrix 

vector of the right hand side of the equations 

work vector 

damping parameters 

stiffness parameters 

test eigenvector 

test measurement location definitions 
numbers of test eigenvectors corresponding to the 
eigenvectors of NASTRAN model 

vector indicating the type of constraints exclusive of simple 
bounds, where IRTYPE(I)=0, 1,2,3 indicates .EQ., .LE., 
.GE„ and range constraints respectively 
vectors containing the lower and upper limits of the general 
constraints 

matrix containing the coefficients of the constraints 



c 

12 

vector containing the coefficients of the objective function 

OBJ 


value of the objective function 

XLB, XUB 

12 

vectors containing the lower and upper bounds on the 
variables 

XSOL, DSOL 

12 

vectors containing the primal and the dual solutions 

ND 


order of the system, default = 63 

N 


order of the system, (input) 

MD 


modes used in the identification, default = 25 

NMODES 


modes used in the identification, (input) 

ID 


number of physical parameters, default =12 

NUNK 


number of unknowns to be identified (input) 

NNK 


number of stiffness unknowns 

NNC 


number of damping unknowns 

NMR 


number of rigid body motion modes 

NMT 


' number of test modes 

LP 


choice of solving techniques 
LP=0, singular value decomposition 
LP=1, constrained optimization 

B. Definition of Input 

and Output Files 

File Name 


Definition 

TEST1, TEST2, TEST5 


files containing 3 different test measurement location 
definitions 

GUY AN_D AMP_F06 


NASTRAN output data file including analytical mass, 
stiffness, damping matrices, eigenvalues and eigenvectors 

GU Y AN_KCOUT 


output file including results 

GUYAN_ELECDAT 


NASTRAN output data file including the grouped element 
matrices Cj 

GUY AN_ELEKD AT 


NASTRAN output data file including the grouped element 
matrices kj 

TEST EIGENV1, .... TEST_EIGENV7 

files containing the test eigenvectors 

TEST_EIGENV AL 


file containing the test eigenvalues 









APPENDIX B 


MSC/NASTRAN INPUT FOR FINAL 
AH-1G SYSTEM IDENTIFICATION RUN 

The following pages include the listings of the input file for the final MSC/NASTRAN runs 
used to compute the structural eigenvalues and eigenvectors for the identified AH-1G structural 
model. 



APPENDIX C 


MSC/NASTRAN OUTPUT FOR FINAL 
AH-1G SYSTEM IDENTIFICATION RUN 

The following pages include the listings of the output files from the final MSC/NASTRAN run 
used to compute the structural eigenvalues and eigenvectors for the identified AH- 1G structural 
model. 




Part II 

IDENTIFICATION OF DAMPING CONSTANTS OTHER THAN THE 
VISCOUS DAMPING CONSTANTS 
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CHAPTER I 
Introduction 

In this part of the work, damping other than viscous damping has been considered. 
As a first step, nonlinear Coulomb damping has been studied. This method can also 
be extended to consider structural damping. This identification procedure uses the 
Hammerstein Feedback Model (HFM) , which represents the nonlinear dynamic sys- 
tem, and Singular Value Decomposition (SVD) Method for estimating the parameters. 
The identification of Coulomb damping constant of a Single Degree of Freedom (SDOF) 
nonlinear dynamic system and the estimating the parameters of Multiple Degree of 
Freedom (MDOF) nonlinear dynamic system have been illustrated in this report. 

The identification of nonlinear dynamical system has received considerably amount 
of attention. These identification procedures are based on various models of nonlinear 
dynamical systems. Usually, a nonlinear system is represented by a set of nonlinear 
differential or integral equations. In many practical applications, an input-output 
approach of a nonlinear dynamical system is a means of describing a relationship 
between the input and the output of the system in some straightforward way and is 
considered to be more useful. 

An approach for modeling a nonlinear dynamical system is by the use of Volt erra Series 

[ 1 ],[ 2 ]. 

x(t ) = [ hi(r)u(t — r)dr 

Jo 

+ / / M T 1 , t 2 M t - T i M < - r 2 )<*Ti dr 2 

Jo Jo 

+ III M T i, T 2,r 3 )tt(< - Tj)u(t - r 2 ) 

Jo Jo Jo 

u(t - T z )dTidT 2 dr 2 + • ■ • 


( 1 . 0 . 1 ) 



The Volterra Series, Eq (1.0.1), is a functional series, It maps past inputs into the 
present output. This means that many kernel values are required to estimate. Several 
techniques have been presented [3], [4], [5]. Because we have to decide which terms of 
Volterra Series are necessary for a given practical problem and to estimate many kernel 
values, the procedure of identification is usually a difficult procedure. 

Several other simple block-oriented input-output, models for representing nonlinear 
dynamical systems are as follows. [7]. 

• Simple Hammerstein Model. 

• Generalized Hammerstein Model. 

• Simple Wiener Approach. 

• Generalized Wiener Approach. 

• Extended Wiener Approach. 

• Simple Wiener-Hammerstein Approach. 

• Generalized Wiener-Hammerstein Model. 

• Extended Wiener-Hammerstein Model. 

The block-oriented models have been widely used because of their simplicity. 

In 1985, a nonlinear difference equation model N ARM AX (Nonlinear Autoregres- 
sive Moving Average Models with inputs) was presented by Leontaritis and Billings 
[9], [10] . The NARMAX model is said to be an unified representation of a finitely 
realizable nonlinear system. The finitely realizable nonlinear system in essence means 
that the state space of the system can not be infinite dimensional. This model maps 
the past inputs and outputs to current output. For the SISO (single input and single 
output) nonlinear dynamical system with white noise, it can be denoted by [11] 

x(k) = F[x(k - 1 ),...,*(* - n x ),u(k - 1), ...,«(* - n„)] (1.0.2) 


2 



Where F(*) is an unknown nonlinear function. In general, it will be determined for a 
given real sampled nonlinear system. Leontaritis and Billings proved that a nonlinear 
discrete time invariant system can always be denoted by Eg. (1.0.2) in a region around 
an equilibrium point, if the response function of system is finitely realizable and a 
linearized model exists at the chosen equilibrium. 

The N ARM AX model is derived assuming zero initial state response, but it can be 
carried over to the non-zero-initial-state cases. The response functions of a system are 
different for different initial condition, but the input-output N ARM AX model for the 
system will always be the same within a region around an equilibrium point. Several 
simple forms of the N ARM AX model have been proposed for nonlinear dynamic system 
identification, such as the Bilinear Model.[ll],[12] 

n n 

x(k) = o 0 -I- ^a,-x(fc - t) + ^2biu(k - i) 

t = l 1 = 1 

+ )«(fc-j) (1.0.3) 

*=i j = i 

the fractional model. [11], [13], [14] 

b[x{k - I),---, 
a[x(k - l),---, 

x(fc - r),u(fc - l),---,t/(fc-r)] (104) 

x(k — r),u(k — 1), • • • ,u(k — r)] 

Haber and Unbehauen [7] prefer the N ARM AX model, because the N ARM AX model 
is parametric and has fewer parameters than the Volterra series. 



In aerospace engineering applications, a nonlinear structural dynamical system is 
usually described by a system of nonlinear differential equations. In SISO case , the 
nonlinear differential equation of a system is of the form 

x -f bx + cx + /(i, x) = u(t) (1.0.5) 
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where /(*) is a nonlinear function of x ,x. If /(*) is represented by a polynomial 
extension for simplicity, Eq.( 1.0.5) becomes 

x + bx + cx + a 2 x + a 3 x + ... + 

&i 2 + fox 3 + ... = u{t) (1.0.6) 

Every term in Eq.(1.0.6) has a distinct physical meaning. Identifying the parameters 
of Eq.( 1.0.6) are useful for dynamic analysis, structural dynamic design, control and 
design modification. If the nonlinear structural dynamic system is modeled by using 
Eq.(1.0.6), the problem of the identification of a system is to estimate the parameters 
. 5, C, Ck 2 , * * ’ ? /^2? " * *• 

Many techniques for estimating these parameters have been proposed. Hanagud, 
Meyyappa and Craig (1985) [15] used the method of multiple scales to formulate a 
procedure for identification of parameters of Eq.( 1.0.6). Mook(1988) [16] used a model 
error method to find the model error d(t) which represents the nonlinear terms of the 
nonlinear dynamic system and then estimated the nonlinear parameters from d(t) by 
using a least square method. Yun and Shinozuka [17] proposed an approach that is 
based on two versions of Kalman filter for identifying the parameters. Ibanez [18] 
used an approach for estimating parameters in which it is assumed that the system 
response is dominated by a periodic response at the forcing frequency and an approx- 
imate transfer function is constructed. Broersen [19] replaced nonlinear terms in the 
equation by a series expansion for a system subjected to random excitation. Distefano 
and Rath, Yun and Shinozuka [20] [21] described several methods of of identification 
and applied nonlinear Kalman filtering techniques for estimation. 

If a structural control is considered, an input- output approach of nonlinear struc- 
tural dynamic system in time domain and its parameter identification is useful. For 
this purpose, the Hammerstein Feedback Model (HFM) has been considered here. 
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CHAPTER II 


Background: Hammerstein Operator and Hammerstein 

Model 


In 1924, P. Uryson investigated a nonlinear integral operator [24], [25] of the following 
type. 

Ax(t) = f k[t,T,x(T)]dT 
Jo 

ten* (2.0.1) 

where Q and SI* are two sets of finite Lebesque measures in a finite dimensional space. 
t 6 SI*, t € S2, — oo < x(t) < oo. K[t,T,x(r )] is measurable and it satisfies the 
Caratheodory condition. The Caratheodory condition is that for all x(t). it is jointly 
measurable in the variables (f,r) € SI* x SI and for all (t,r) € S2* x SI, it is continuous 
in :e(t). 

In 1930, A. Hammerstein studied the following integral equation: 

x(t)= f K(t,T)f(T,x(T)]dT (2.0.2) 

Jo 

This kind of equation is known as Hammerstein equation. 

Eq.(2.0.2) is a special form of the Uryson nonlinear integral operator : 

/ Ko(t,T)f[T,x(T)]dr = Hx(t) (2.0.3) 

Jo 

This integral operator H is called Hammerstein integral operator. Hammerstein inte- 
gral operator H can be denoted by the following form: 


H = Ko? 


(2.0.4) 



where A'o represents a linear integral operator with kernel Ko{t,r): 

K 0 x(t) = [ Ko(t, t )x(r )dr (2.0.5) 

Jn 

and T represents the nonlinear superposition operator [25]. 

Tx{t) = /[t,x(t)] (2.0.6) 

Then Hammerstein equation Eq.(2.0.2) can be expressed by 

x(t) = K<>Tx(t) (2.0.7) 

and Hammerstein integral operator can be denoted in following form 

Hx(t) = K Q Fx{t) (2.0.8) 


Let L q ,X ( 3,X7 express the sets of measurable function x(r) . They have the norms 
separately as follows. 


II X IU= 1 

:/ 1 

J 0 

(2.0.9) 

II 1 b= \ 

[/ \ x(r) r e drf 
J n 

(2.0.10) 

II 1 Ilf = 

and Lo is the set of z(r) with norm 

\j n 1 *(t) r *-r 

(2.0.11) 

II * ll« 

3 = sup I x(r) I 

(2.0.12) 


There are two theorems about Hammerstein operator H: [25] 

Theorem 1: Let T be an operator acting from L a to L (7 > 0) and let A 0 be a 
continuous operator acting from L y to L$. Then Hammerstein operator H = ho^F 
acts from L a to Lp and is continuous. 

Theorem 2: Let T act from L a to L 0 and let A' 0 be a regular operator acting from L 0 
to Lp (0 > 0), then the Hammerstein operator H = KoJ- acts from L a to L$ and is 
continuous. 
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The theorems 1,2 are very useful. It permit one to construct the Hammerstein operator 
if A'o and T are known. 

In 1966, K.S. Narendra and P.G. Gallman [6] suggested a Hammerstein Model for 
identification of nonlinear dynamic system. They assumed that the response x{t) of 
nonlinear dynamical system is 


x(t) = Hu(t) 

= K 0 Tu{t) (2.0.13) 

where u(t) is an input function. Then the Hammerstein Model suggested by Naren- 
dra and Gallman consists of a nonmemory (independent from history) nonlinear gain 
having a polynomial form followed b}- a linear discrete system. The nonlinear gain is 

v(t) = Tu(t) 

= ciu + c 2 u 2 + 1 - Cj,u p (2.0.14) 


The linear discrete system has a response . 


x(t) = / Ko(t,T)v(r )dr 
J n 

Narendra and Gallman used the following form to denote the linear system. 

biq 1 -f • ■ • + b n q 
a 0 + aig -1 H 1- a n q~ n 

where the q' 1 , • • • ,g - " are the time delay operators. They are defined as 


(2.0.15) 


(2.0.16) 


x(k)q 1 = x(k — 1) 


x(k)q n = x(k — n) (2.0.17) 

This Hammerstein Model suggested by Narendra and Gallman is illustrated in Fig. 2.1. 
Hammerstein model provides a simple input-output model for identification of nonlin- 
ear dynamical system. In the past years , Hammerstein model has been widely used 
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Figure 2.1: The Hammerstein Model 


in various fields. 


Neil, Francis and Rein [26] presented a simple iterative technique for estimating param- 
eters in a Hammerstein model for a case where noise in the output data is correlated. 
They suggested a Hammerstein model with noise . The noise is modeled by the fol- 


lowing transfer form . 


H(q) = 


(2.0.18) 


l + c 1 f 1 +-- + c,g- ra 

The expected value of noise d{j) equals to zero. The parameters are estimated by an 


iterative method. 

Greblick and Pawlak [27] used a correlation method for Hammerstein system identifi- 
cation by using non-parametric regression estimation. Shih and Kung [30j, [31 j used 
shifted Legendre and Chebvshev expansions for identification of the nonlinear dynamic 
system described by a Hammerstein model. The input and output are expressed b\ 
using a finite number of the shifted Legendre polynomials or of the shifted Chebyshet 
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polynomials. Substituting these into the Hammerstein model in a state variable form, 
and integrating from 0 to t, then we reduce the nonlinear equations to a linear alge- 
braic equation. The parameters are obtained from the linear algebraic equations. 
Horng and Chou, [32], used shifted Jacobi series to express the input and output. Sub- 
stituting the input and output into the Hammerstein model, then integrating the 
model from 0 to t, a linear algebraic equation is obtained for identification of nonlin- 
ear dynamical system. Chung and Sun, [33], used a Taylor’s series approximation for 
Hammerstein model to estimate the parameters. Kung and Shih [34], and Jiang [35], 
used the Block pulse function for identification of parameters of a nonlinear system 
with Hammerstein model etc. 

Actually, the Hammerstein model suggested by Narendra and Gallman can be consid- 
ered as a superposition of several linear system models. Their objective was primarilly 
to consider only nonlinearity in the forcing function u(f), however in a structural dy- 
namic system, nonlinearities are from damping and stiffness terms, (or the plant model 
). Therefore we will propose a feedback model which is named here as Hammerstein 
Feedback Model (HFM) for identification of nonlinear feedback system for a nonlinear 
structural dynamical system and identify the nonlinear plant parameters. 
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CHAPTER III 

HAMMERSTEIN FEEDBACK MODEL OF NONLINEAR 

DYNAMIC SYSTEM 


Consider a linear equation. 

Lx = /(f) (3.0.1) 

where f(t) is an almost periodic function (ap-function). The ap-function is defined as: 
A function x(t) E C(R) ( C(R ) denotes a continuous metric function space ) , which 
has the translate 

Xh{t) = x(t + h) (3.0.2) 

where ( — oc < i < oo) and ( — oo < h < oo), is called as an ap-function, if its translates 
form a compact set in C(R). 

L is a regular ap-operator . 

L = IF + + ” ' + An ^ ^ 3 '°‘ 3) 

where Ai(t) (i = 1, •••,«) are ap-functions. There is at most one function G(t,r) 
( — cxd < t,r < oo) , which is called to be the Green’s function, such that one can write 
the solution of Eq.( 3.0.1) as: 

x(<)= / G(t,r)/(r)<£r (3.0.4) 

Jn 

Now consider the nonlinear equation. 

Lx := f(t,x,x,u) (3.0.5) 

where the function f(t,x,x,u ) is jointly continuous for t, x and almost periodic in t . 

It can be easily seen that the ap-function x{t) is a solution if and only if it is a solution 



of the integral equation, called the Hammerstein integral equation. 


x{t) = [ G(t,r)/[r, x(r),i(r),u(r)]<fr 

Jo 

= H(x,x,u,t ) (3.0.6) 

The right side of Eq.(3.0.6) is known as Hammerstein integral operator, which can he 
denoted by 

H = K 0 T (3.0.7) 

where Kq represents a linear integral operator with Green’s function G of operator L, 
and T represents a nonlinear superposition operator 

F{x,x,u,t) = f[x{t),x{t),u{t),t ] (3.0.8) 


Derivation of Hammerstein integral equation 

For linear case, consider a linear dynamic system for instant. 

x(t) = Ax(t) + Bu(t ) 

where A and B are constants. The Laplace transform of Eq.(3.0.9) is 

sx + x(0) = Ax + Bu 


i.e. 

x = [si — j4] -1 x(0) + [si — A}~^Bu 
The response of system is 

x(t) = £^^(0) + [' e At e~ Ar Bu{T)dT 
Jo 

i.e. 

x(t) = [ GBu(t)cIt 
Jo 


(3.0.9) 


(3.0.10) 


(3.0.11) 


(3.0.12) 


(3.0.13) 
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If the initial condition is assumed as x(0) = 0. For nonlinear case, consider a nonlinear 
dynamic system 

i(t) = Ax(t) + /[x(f),i(f ),«(<)] (3.0.14) 

Initial condition: x(0) = 0. Construct a response for the nonlinear dynamic system, 
which is denoted by Eq.(3.0.14). 

x(t) = [ d At e~ Ar /[x(r),x(r),u(T)]dr (3.0.15) 

Jo 

The derivative of Eq.(3.0.15) is 

x{t) = J o e~ AT flx{T),x{T),u(T))dr 

+ e At e~ At f[x(t ), i(<), u(t)] = Ax(t) ■+ f[x, i, u] (3.0.16) 

Eq.(3.0.15) is the Hammerstein integral equation, which is equivalent to the nonlinear 
differential equation Eq.(3.0.14). Eq.(3.0.15) can be rewritten as 

x{t) = H(z,i,u) (3.0.17) 

where H is the Hammerstein integral operator. Since there are two theorems about 
Hammerstein operator H : 

Theorem 1 : Let T be an operator acting from space L Q to space L ^ (7 > 0) and let 
A"o be a continuous operator acting from space L 7 to space Lp. Then Hammerstein 
operator H = KoJ~ acts from space L Q to space Lp and is continuous. 

Theorem 2: Let T act from space L a to space L 0 and let Ko be a regular operator 
acting from space Lo to space L$ (/ 3 > 0), then the Hammerstein operator H = K$T 
acts from space L a to space Lp and is continuous. 

According to Hammerstein integral equation and theorems 1 and 2 of Hammer- 
stein integral operator, a Hammerstein Feedback Model (HFM) can be constructed for 
representing a nonlinear dynamic system. The HFM consists of nonlinear part, ’which 
is expressed by a superposition operator T and contains the nonlinear terms of state 
variables, followed by a linear system, which is denoted by linear integral operator 
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Figure 3.1: The Hammerstein Feedback Model (HFM) 

A' 0 . HFM is a simple block-oriented input-output model for identification of nonlinear 
structural dynamic systems. It can be illustrated in Fig. 3.1. 
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CHAPTER IV 

IDENTIFICATION OF NONLINEAR STRUCTURAL 
SINGLE DEGREE OF FREEDOM (SDOF) DYNAMIC 
SYSTEM BY USING HFM 


In chapter 3, the Hammerstein Feedback Model (HFM) of nonlinear dynamic sys- 
tem has been proposed for identification of nonlinear dynamic systems. In this chap- 
ter the HFM is used for identification of nonlinear structural single degree of freedom 
(SDOF) dynamic system. The HFM in discrete time domain of SDOF nonlinear struc- 
tural dynamic system has been derived. The Singular Value Decomposition Method 
(SVDM) is used for estimating parameters of HFM of SDOF nonlinear structural dy- 
namic system. 


4.1 Hammerstein Feedback Model in Discrete Time Domain of SDOF 
Nonlinear Dynamic System 

For HFM, the response of nonlinear dynamic system is the convolution of a weighting 
function i.e. Green’s function of linear dynamical system and a nonlinear function of 
the input and the output of nonlinear dynamic system. The Z-transformation of a 
sequence of function and its properties can be used for deriving the HFM in discrete 
time domain. 


Consider a sequence of function f(k) ( f(k ) = 0 for k = — 1, — 2, • • •)- The Z- 



transform of f(k) is defined by 


F(-) = E /(*)•'* <4-1.1) 

k = 0 

and inverse Z-transform is 

f(k) = Z-'\F(.z)} (4.1.2) 

where z is an arbitrary complex number. The Z-transformation of a sequence of func- 
tion f(k) has following properties: 

Linearty: 

If f{k) = afi(k) + bf 2 (k) for k = 0,1,2, •••, where a, b are constants, then 

F(z) = aF 1 (z) + bF 2 (z) (4.1.3) 

for | z |> max(/?!, J? 2 ), where R u R 2 denote radii of convergence for Fi(z), F 2 (z), 
respectively. 

Right-shifting Property: 

Consider f{k) ( f{k ) = 0, for k = 0, -1, -2, • • •) and y(k) = f(k-m) for k = 0, 1, 2. - - 
From the definition of Z-transformation, we have 

y(r) = f^f(k-m)z~ k 

k-0 

= + 1 - m)z~ 1 + • • • + /(-l)z- m+1 

- f/(0)z- m + /(l)r- m+1 + ••• 

= r- m [/(0) + /(l)z- 1 + ---j 
= =- m F(z) 

for \z\> R (4.1.4) 

where R is the radius of convergence for the Z-transformation of f{k). Then the 
Right -shifting property of Z-transformation can be denoted by the following form. 


Z[f(k-m)) = z- m F(z) 


(4.1.5) 



Convolution-Summation Property: 

Consider a convolution y(k) of two sequences h(k) and u(k) 


y(k) = ~ *) 

»=0 


The Z-transformation of y(k) is 


y ( : ) = EE fc (» M fc -*)] 2 k 

k = 0 »= 0 

From the Right-shifting property, we have 


Z[u{k-i)\ = z~ i U{z) 


(4.1.6) 


(4.1.7) 


(4.1.8) 


Eq.(4.1.7) can be represented as 


Y(z) 




H(z)U(z) 


(4.1.9) 


If input sequence u(k ) is kroneker 6 function, 


u(k) = ^(fc) 


(4.1.10) 


the Z-transformation of u(k) is 


U(z) = 1 


(4.1.11) 


In this case, Eq.(4.1.9) becomes 

Y(z) = H(z) (4.1.12) 

Eq.(4.1.12) denotes that response of any linear discrete system for Kronecker 6 input 
is equal to the weighting sequence of the linear dynamic system. 


The Hammerstein Feedback Model implys that a nonlinear dynamic system is as- 
sumed as a linear system with a nonlinear input, which is a function of the input 
and the output of nonlinear dynamic system. For a nonlinear dynamic system, if the 


16 



responses of system are known and steady, and if the parameters of system are time- 
invariant and linear, the estimating parameters of system reduces to a linear problem. 


Let us consider a SDOF linear dynamic system. This SDOF linear dynamic system 
is represented in a differential equation with order n in general. 


d"x(t) 

dt" 


+ j#„_i + • ' ' + 


= &i 


dt' 

d n ~ 1 u(t) 

dt”" 1 


+ 1- b„u(t ) 


(4.1.13) 


with initial conditions: x(0), i(0), •••, If the initial conditions are assumed 

to be zero, 


z(0) = 0 


i(0) = 0 


d" _1 

ric(O) = 0 

dt" -1 

the response of linear dynamic system has the response of following form 

z(t) = [ h(t — r)u(r)dr 
J n 

In discrete time domain, the sequence of response of system is 

x(k) = ^2 h(k)u(k — i) 

1=0 

The transfer function of linear system in z-domain is 

X{z) 


H{z) = 


U(z) 

f>i z" 1 4- ... 4- b n 


z" + aiz" -1 -f ... -t- a„ 

The transfer function Eq.(4.1.17) can be rewritten as 

bj z 1 + • • • 4- b„ z 


(4.1.14) 


(4.1.15) 


(4.1.16) 


H(z) = 


1 + ajz 1 + 1 - a„: 


(4.1.17) 


(4.1.18) 
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Consider a nonlinear superposition operator J~ > It is assumed in the form of a 


polynomial expression for simplicity. 

T [ z(f),x(f),ti(f),i] = + 72T 2 (0 + • • • 

+ 7 p r p {t) + tt 2 x 2 (t) + (4.1.19) 

Since the function f(x,x,u,t), which is respect with the nonlinear superposition op- 
erator, is a polynomial expression, it satisfies the Caratheodory condition. The z- 
transform of Eq.(4.1.19) is 

Z{y) = 'fiZ(v) + 72Z(x 2 ) + • • - + 7 P Z(x p ) 

+ p 2 Z(i 2 )-f ... + /*,Z(x*) (4.1.20) 


According to Eq.(4.1.9), we can construct a Haromerstein Feedback Model in dicrete 
time domain, which is 


Z(x) 


+ 


hz- 1 + ••• + 

1 + OiZ -1 + •••-+- CL n z~ n 


[7i Z{u) 


72 Z(x 2 ) + • • * + 7 P Z(x p ) 


+ p 2 Z(i 2 ) + + fi q Z(x q )] 


(4.1.21) 


This HFM of SDOF nonlinear dynamic system in Z-domain is illustrated in Fig.4.1. 
Eq.(4.1.21) can be rewritten as 

Z(x) — —a- l Z(x)z~ 1 — --- — a„Z(x)z n + b 1 j 1 Z(v)z 

+ 6 1 7 ,Z(* a )z - 1 + ••• + b llp Z(x p )z~ l + b^ 2 Z(i 7 )z^ 

+ • • • + b- l fi q Z(x q )z~ 1 + h b n iiZ{u)z~ n 

+ ■■■ + b nf i q Z(x*)z-'' (4-1.22) 


By using the Right-shifting Property of Z-transform, Eq.(4.1.22) yields 
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Figure 4.1: The HFM of SDOF nonlinear dynamic system in z-domain. 

Z[x(k)} = —a\Z[x{k — 1)] o n Z[*(fc-n)] + 6i7i^Mfc- 1)] 

4 br^Z[x 2 (k - 1)] 4 • • ■ + b plp Z[x p {k - 1)] 4 b lt i 2 Z\x\k - 1)] 

4 • • • 4 - 1 )] + ••• + b n ^Z[u{k - n)] 

4 • • ■ 4 b n fi q Z[x q (k - n)] (4.1.23) 

Consider the inverse Z-transform, we have 

i(fc) = — aix(k — 1) — a„x(k — n) 

4 6ni«(fc - 1) + fci723- 2 (A- - 1) 

-f • • • 4 bi~) p x p (k — 1) 4 bifi2X 2 (k — 1 ) 

-I- b b„-)Mk - ri) + b„') 7 T 2 (k - n) 

4 h b n ^ p x p (k - tj) 4 b„n 7 i i {k ~ * ) 

4 rb n fi q x{k-n) (4.1.24) 

Eq.(4.1.24) is the difference form of HFM in discrete time domain of SDOF nonlinear 
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dynamic system. This HFM in discrete time domain, Eq.(4.1.24), is equivalent to 
following nonlinear differential equation. 


d n x(t) 

dt n 


cT 1 x{t) . . 

+ fli — 77—; h 1- a n z(0 


= b 


dt n ~ x 
dt n ~ l 


+ ”- + 6 n /(0 


where input f(t) is 


/(<) = 7i«(0 + 72^ 2 (0 + "- 

+ Kp2 P (t) + V2X 2 (t)+ + Hqi q (t) 


(4.1.25) 


(4.1.26) 


Observably, since the parameters of nonlinear differential equation, Eq.(4.1.25), are 
indepedent from the initial conditions, the parameters of HFM in dicrete time domain, 
Eq.(4.1.24), are independent from the initial conditions. This means that we can use 
the input and output data under any initial conditions to identify the parameters of 
HFM in discrete time domain of a nonlinear dynamic system. 


4.2 The HFM of a Nonlinear Structural SDOF Dynamical System 

A nonlinear structural SDOF dynamical system usually is expressed in the form of 

x + bx *f cx + f(x, x) = u (4.2.27) 

where f(*) is a nonlinear function of x,x. It can be also approximately expressed by 

x + bx + cx + Q 22: 2 + <* 3 £ 3 + • • - 
+a p x p + /3 2 i 2 + f 0 q x 9 = u (4.2.28) 

In this case, the linear differential equation corresponding to Eq.(4.2.27) or Eq.(4.2.28) 
is 

x + bx + cx = u (4.2.29) 
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with initial condition: x(0),x(0) 

The transfer function of the linear system in z-domain is 

X(z) 


H(z) = 


U(z) 


.-2 


1 4 a,\z 1 4 fl- 2 z 2 

The nonlinear input f(t) in HFM is defined as 

f(t) = u(t) 4 7 2 x 2 (t) 4 f7 p xP (t) 


4 P2^ 2 4* • • ■ + 


(4.2.30) 


(4.2.31) 


The z-transform of nonlinear input f(t) is 

Z{f) = z{u) + 'r 2 Z{x 2 ) + --- + i P z{x p ) 

4 p 2 Z{x 2 ) + *r fi q Z{x q ) (4.2.32) 

Then we have the HFM of nonlinear structural SDOF dynamic system in z-domain. 


— 2 


Z(x) = 


— 2 {Z{u) 


1 4 oiz -1 4 a 2 - 
4 'y 2 Z(x 2 ) 4 • • • + lpZ{x p )} 

4 y. 2 Z{x 2 ) 4 • • • + PqZ(x v )} 


(4.2.33) 


Eq.(4.2.33) can be rewritten as 

Z(x) = -aiZ(x)z~ 1 - a 2 Z(x)z~ 2 4 Z(u)z~ 2 
4 72 Z(x 2 )z~ 2 + (-7 pZ(x p )z 

4 fi 2 Z(x 2 )z~ 2 4 1 - fi q Z{x q )z~ 2 (4.2.34) 

By using the Right-shifting Property of Z transform , we have the HFM of nonlinear 
structural SDOF dynamic system in discrete time domain. 

x(k) = — a\x(k — 1) — a 2 x(k — 2) 4 u(k — 2) 

4 a 3 x 2 (k- 2)4 f a p+l x p {k - 2) 

4 a p + 2 x 2 (k — 2) 4 • • • + a p +g+-[i q (k — 2) (4.2.35) 
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where 


a 3 = 72 


^4 = 73 


a p+q 41 — Vq (4.2.36) 

The relationship between parameters of nonlinear differential equation and param- 
eters of HFM are 

b — (di -f- 2)/At 

c = (oi + a 2 + 1)/(A f) 2 

a 2 = -a 3 /(Af) 2 

q 3 — —a^KAt) 2 

02 = — a P +2/(Af) 2 

Pq = - Wl /(A*) 2 (4-2.37) 

4.3 Estimation of Parameters of Nonlinear Structural SDOF Dynamic Sys- 
tem 

When a nonlinear dynamical system is modeled by HFM in discrete time domain, 
Eq.(4.1.24), and N + n samples of the input and the steady output from k — n to k + N 
are substituted into Eq.(4.1.24), we have N equations. The set of equations can be 
represented in matrix as the following form. 

X = >10 (4.3.38) 
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where 


*(*) 

x(k + 1) 

x(k + N) 


(4.3.39) 


— x(k — 1) 

£ 

1 

H 

1 

u(k — 1) 

-x(k) 

• ■ • — x(k — n -f 1) 

u(k) 

— x(k + A 7 - 1) 

• • • — x(k -f N — n) 

u(k + N 

■ • x 2 (k — 1) 

• • • x p (k — n) 



••• x\k ) ••• z p (fc-n + l) 


••• x*(k + N- 1) ••• x p {k + N-n) 

x 2 (k — n) ••• x g (k — n) 

x 2 (k — n + 1) x 9 (k — n + 1) 

x 2 {k + N-n) ••• i g (k + N-n). 


(4.3.40) 


ai 

an 

£>i7i 


(4.3.41) 


Solving the Eq.(4.3.38) we will obtain the parameter vector 0, which dominate the 
nonlinear dynamical system. 

For identification of nonlinear structural dynamic system, the structural HFM in dis- 
crete time domain is considered. When N + n samples of input and output are taken, 
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the Eq.(4.2.35) becomes the following form in matrix. 



II 


— x(k — 1) 

-x{k- 2) 

u(k — 2) 

”*(*) 

—x(k — 1) 

u(k — 1) 

— x(k + N 

-1) -x{k + N 

— 2) t /( k -+• N 


x 2 (k — 2) 

■■■ x p (k- 2) 


x 2 (k — 1) 

••• x p (k — 1) 


x 2 {k + N -2) 

x p (k + N 


i 2 {k- 2) 

••• x q (k- 2) 


x 2 {k- 1) 

• • • x 9 (k — 1) 


x { k AN -2) 

• •• x 9 (k + N 


(4.3.42) 


(4.3.43) 


0 = 


a i 
a 2 

O-p - 1 


(4.3.44) 


a p+q + 1 

A modeling error noise e(fc) usually is considered. This noise is assumed to be white. 
In this case, Eq.(4.3.38) is 

X = A© + e (4.3.45) 


Estimating the parameter vector © from Eq. (4.3.45) is a standard least squares prob- 
lem. The problem is to identify (estimate) the parameter vector © which minimizes 
the || A@ - X || . There are several methods available for estimating parameters from 
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Eq.(4.3.45). The basic method for estimating © is the least squares method. The 
parameter vector can be estimated by 

© = (A t A)~ 1 A t X (4.3.46) 

Since matrix A of HFM has elements x‘(fc - j) (i - l,2,...,p; j = l,2,...,n), the 
determinant of some submatrix of matrix A is the Vandermonde determinant. 

x 2 {k — 2) x p (k — 2) 

x 2 (k — 1) ••• x p {k- 1) 

det . 

r 2 (Hp- 2) ••• x p (k + p — 2) 

= x 2 (k — 2) ■ ■ • x 2 (k + p — 2) 

1 ... x p ~ 2 (k- 2) 

1 ... x p ~ 2 (k- 1) 

1 ••• x p - 2 (k + p- 2) 

The Vandermonde determinant has following value. 

= n (4 - 3 - 48) 

l<j<t<n 

if At is a very short time period, the difference between x{k - j) and x(k - j + 1) is 
small value, then the value of determinant of maximum square submatrix of matrix A 
is small. The matrix A becomes an ill-conditioned matrix. Any small round-off errors 
in the elements of A can cause large changes in {A T A) 1 . In this case, Eq.(4.3.46) is 
not stable for identification of HFM. 




(4.3.47) 
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4.4 Singular Value Decomposition Method 


In this report, the Singular Value Decomposition (SVD) method has been used for 
estimating the parameters of HFM. For an arbitrary matrix A (A € R mxn ), if there 
a matrix B which lets AB , BA be a H matrix ( H matrix: m x n , m > n, and 
rank(H) — n.) and 

ABA = A (4.4.49) 

BAB = B (4.4.50) 

the matrix B is defined as a pseudo-inverse matrix A + of A . 

There is a theorem of singular value decomposition (SVD) for an arbitrary matrix 

A. 

Theorem: if A 6 ? there exist orthogonal matrices 

U — [u i , • • ■ , ix m ] 

IT £ R™*™ (4.4.51) 

and 

V = [Vi, •••,!>„] 

V € ir* n (4.4.52) 

such that 

U t AV = diag(s u - ••,»,) (4.4.53) 

where Si > ^2 > ■ • • > s p > 0. According to this theorem, any matrix A (m x n) with 
rank r can be decomposed to 

A = USV T (4.4.54) 

where U(m x m), V(n x n) are orthogonal matrices. S is a m x n matrix with all 
elements of which is equal to zero with the exception of the first r diagonal elements 

Si > $2 > • • • > s T > 0 (4.4.55) 
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We can construct a pseudo-inverse matrix A + of A as 


A + = VS + U T 


(4.4.56) 


where S + is a n x m matrix and has all zero elements except the first r diagonal 
elements. The nonzero elements of S + are ’ s ^ 1- This A + satisfies the 

definition of pseudo-inverse matrix. 

A standard least squares problem is finding a vector € R n for equation >10 = A 
where A € R m * n and A" € R m and m > n, such that 


min || >10 — X H 2 (4.4.57) 

Here || >10 — A” H 2 is the p-norm ( p = 2) of vector. It is defined as 

|| A0 - A' || 2 = [(A© - A') r (A0 - A')]* (4.4.58) 

Denote the minimum sum of squares by p\ t 

p], =|| AQu ~ X \\l (4.4.59) 

When a matrix >1 (m x n) is decomposed by using the orthogonal matrices V and F, 
we have 

||A0-A'||’ = \\U T AV(V T Q)i-U T X \\l 

= j;Mv r e),-«*W 
+ x>w 

i= r+1 

Observably, if © has the following form. 

@ lt = £(vJX/s,)v r 

i = 1 

then the minimum of || ^4© — -V ||| is 

A = £ («W 

i = r + 1 


(4.4.60) 


(4.4.61) 


(4.4.62) 
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From Eq. (4.4.61) and Eq. (4.4.56), we obtain the solution of least squares problem. 

= A + X (4.4.63) 

where A + = VS + U T . Actually, if the rank of matrix A, r is equal to n, then A + can 
be (A t A)~ 1 A t and the least squares solution 0j, is 

Q h = {A t A)~ 1 A t X (4.4.64) 

In this report, the U, V and singular values Si of matrix A are calculated by using the 
program of IMSL. 


4.5 Numerical Examples 

Several numerical examples have been presented by using HFM for identification of 
nonlinear structural SDOF dynamic system in this section. The nonlinear structural 
dynamic system are described by nonlinear differential equation. The response x(t ) and 
velocity x(t) are obtained by using Runge-Kutta Method. The estimated parameters 
of nonlinear terms are compared with the real parameters. 


Example 1: 

The simplest case has been considered in this example. A nonlinear dynamic system 
is assumed to have the following differential equation. 

2.56a + 0.32a + x + 0.05a; 3 = 2.bcost (4.5.65) 

It is from Mook’s paper. In Mook’s paper, [16] a linear model is assumed as 

2.56a + 0.32a + x = 2.bcost (4.5.66) 

The real nonlinear dynamic system can be represented as 

2.56a + 0.32a + a = 2.bcost + d(<) (4.5.67) 
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where d(t) is model error. The model error d(t) between nonlinear dynamic system 
and assumed linear model is estimated by using Two Point Boundary Value Problem 
(TPBVP) method. Then the model error is assumed to consist of two nonlinear terms. 

d{t) = ax 2 (t) + 0x\t) (4.5.68) 

The parameters a and 0 are estimated by the least squares method from the model 
error d(k) in discrete time domain. 

Since this is a nonlinear dynamic system with one order, the linear dynamic system 
considered for constructing HFM has following transfer function in z-domain. 

H(z) = ■ (4.5.69) 

1 -I- a\z 1 

The nonlinear input of HFM is assumed as 

f{t) = 2.5 cos(t) + 72 * 2 (<) + 73 * 3 (t) (4.5.70) 

Then we can assume a HFM for the nonlinear dynamic system. 

x(k) = —aix(k — 1) + a 2 x 2 (k — 1) -f a 3 x 3 (k — 1) (4.5.71) 

The response of system is obtained by Runge-Kutta method from Eq.(4.5.65). The 
input and response are considered as data for identification of this nonlinear dynamic 
system by using HFM. 

Sampling period At is assumed to be 0.01. There are 628 samples taken in a period. 
901 samples of input and output of system are used for forming the equation of least 
squares problem. The parameters a*, a 2 , a 3 are estimated by using SVD method. The 
nonlinear estimated parameters a, (3 are obtained from equations: 

a = a 2 /(A<) 2 

0 = a 3 /(At) 2 (4.5.72) 

They are compared with mook’s results in table 4.1. 


29 



Real P. 

Est. P. by Mook 

Est. P. by HFM 

0 

0.00001 

0.00002 

0.05 

0.0492(error:1.6%) 

0.04995( error :0.09%) 


Table 4.1: The estimated nonlinear parameters 


Example 2: 

A nonlinear SDOF structural dynamic system described by well-known Doffing’s 
equation is considered for this example. The Duffing's equation has the following form. 

x + 0.225x 4- 0.0025x 3 = OMcosO.bt (4.5.73) 


The initial conditions are assumed as 


x(0) = 4.0 

x(0) = 0 (4.5.74) 


The response x(t) of Eq.( 4.5.73) is calculated by using the Runge-Kutta method and 
illustrated in Fig. 4.2. 

The linear dynamic system considered fo HFM has the transfer function in z- 
domain. 


H(z) 


1 -|- a 2 z 1 + a 2 z 2 


The nonlinear input of HFM is assumed as 


(4.5.75) 


f(t) = 0.02cos0.5(0 4- i 2 X 2 (t) - 7s x 3 (t) (4.5.76) 


The HFM of this system is assumed as 

x{k) - -a^k - 1) - a 3 x(k - 2) - a 3 x 7 (k - 2) 
-f a A x 3 {k — 2) 4- 0.02cos0.5tA- - 2) 


902 samples of input and output of system are used for parameter estimation. The 
sampling time period A / is assumed as A / = 0.05. The parameters of HFM are 
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displacement 



time (sec.) 


Figure 4.2: The response of a SDOF nonlinear dynamic system 
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Exact P. 

Estimated P. 

Error 

0.0 

0.000012 


0.0025 

0.0025030 

0.12% 


Table 4.2: The estimated nonlinear parameters 

estimated by using SVD method. 

ai = -1.999943 

a2 = 0.999995 

a 3 = 0.320225 x 10" 6 

a 4 = -0.625760 x 10 -5 (4.5.78) 

The parameters of nonlinear terms of nonlinear differential equation are obtained by- 
using Eq.(4.2.37) and shown in table 4.2. 

Example 3: The third example is the Duffing’s equation, which has linear damping 
term. This example can be expressed as 

x + O.Oli + 0.225i + 0.0025* 3 = 0.02cos0.5t (4.5.79) 

The initial conditions are 

x(0) = 4.0 

i(0) = 0 (4.5.80) 

The response of this system is illustrated in Fig. 4.3. 

The HFM is the same with in Example 2. That is 

X{k) = -a 1 x{k-l)-a 2 x(k-2) + a 3 x 2 (k-2) 

+ a 4 x 3 {k - 2) + 0.02cos0.5(fc - 2) (4.5.81) 

In the HFM of this example, a quadric nonlinear term and cubic nonlinear terms are 
assumed. 
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displacement 






Exact P. 

Estimated P. 

Error 

0.0 

0.0000064 


0.0025 

0.00256 

2.4% 


Table 4.3: The estimated nonlinear parameters 

Sampling time period At = 0.05 is considered and 755 samples have been used. 
Estimated sample range approximately is three period of the input and the output. 
The results obtained are 


a a = -1.99890 

a 2 = 0.999460 

a 3 = 0.162931 x 10 -7 

a 4 = -0.6407975 x 10" 5 (4.5.82) 

The estimated parameters of nonlinear terms of Duffing’s equation, Eq.(4.5.79). are 
listed in table 4.3. The results show that there is not the quadric nonlinearity of 
displacement x(t). 


Example 4: 

Coulomb damping is the friction force between two contact surface. There is a 
friction drag force. 

c | x | x (4.5.83) 

where c is a positive constant. A nonlinear SDOF structural dynamic system with 
Coulomb damping is considered. This nonlinear dynamic system is denoted as following 
differential equation. 

x + 0.225x + 0.1 | x | x = 0.02 cos 0.5f (4.5.84) 
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The initial conditions are 


x(0) = 4.0 

i(0) = 0 


(4.5.85) 


The responses x{t), x(t) of Eq.(4.5.84) are calculated by using Runge-Kutta method 
and illustrated in Fig. 4. 4. Sampling time period At is assumed to be 0.05 and 755 
samples of input and output are considered as data for identification. 

The HFM of this system is assumed as 


x(k) = -aix{k - 1) - a 2 x(k - 2) + a 3 | i(k - 2) | x(k - 2) (4.5.86) 


Then a least squares problem is formed for estimating the parameters < 21 , 62 , < 23 * The 
damping constant c of Eq.(4.5.84) can be obtained by 


c = 


03 

At 2 


(4.5.87) 


The parameters, a j, 02 , 03 of HFM are estimated by SVD method They are 


a, = -1.99937 
a 2 = 0.999926 
a 3 = -0.00024675 

The Coulomb damping constant c estimated is 0.0987004. 


(4.5.88) 

The error is 1.29%. 


Example 5: 

Usually a linear SDOF structural dynamic system 

mx + cx -f kx = u (4.5.89) 

can be modeled by a discrete difference model for identification of system. (Fig. 4. 5) 
The difference model is 

x(fc) = -a-[x{k — 1) — Q 2 x(k — 2) + u(k) (4.5.90) 
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Figure 4.5: The linear dynamic system 

where oj and a 2 are unknown parameters, which dominate the linear dynamic system. 
The Model, Eq.(4.5.90) , is an input-output approach of linear dynamic system. It is 
not possible to obtain the parameters m, fc, c in Eq.(4.5.89) from the parameters a-j, 
a 2 . 

In engineering, the m, k are easy to be obtained from the real structure. The damping 
constant c is necessary to estimate. In this case, the Hammerstein Feedback Model can 
be used for identification of damping constant c. The Unear dynamic system is assumed 
as a feedback Unear dynamic system. It is iUustrated in Fig. 4. 6. The parameters Oi 
and a 2 are calculated from m. k. The HFM of the system is to be 

x(k) 4 a\x(k - 1) - a 7 x(k - 2) = a 3 x{k - 2) - v{k - 2) (4.5.91) 

Consider a Unear dynamic system 

x + 0.225.r * 0.0 lx = 0.02 cos 0.5/ (4.5.92) 
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Figure 4.6: The linear model with damping feedback 


with initial conditions: 

x(0) = 1.0 

i(0) = 0 (4.5.93) 

The responses x(f) and i(t) are illustrated in Fig.4.7. The sampling time period 
Al = 0.05 and Too samples of input and output are taken for identification. The 
estimated damping constant c is 0.00999692. The error is 0.03 /o. 
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Figure 4.7: The responses 
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CHAPTER V 

IDENTIFICATION OF NONLINEAR MDOF 
DYNAMICAL SYSTEM BY USING HFM 


A large and complex structural system is usually approximated by a multiple de- 
grees of freedom (MDOF) dynamical system by using method, like the finite element 
method. The identification and modeling of nonlinear MDOF dynamical systems by 
the use of input and output data is then a very important problem in practical struc- 
tural dynamical system. Masri, Millar, Sand, and Caughey [42] [43] presented a self- 
starting multistage time-domain procedure for the identification of nonlinear MDOF 
dynamical systems in free oscillations or subjected to an arbitrary direct force ex- 
citations and nonuniform support motions. Yasuda, Kawamura and Watanabe [22] 
presented a technique in frequency domain for identification of nonlinear MDOF dy- 
namical system. This technique is as follows. The periodic steady state responses data 
are measured from a MDOF nonlinear dynamical system subjected to a periodic force 
excitation. The nonlinear terms are expressed in terms of polynomials with unknown 
coefficients. The parameters are determined by expressing the quantities in a Fourier 
series and by applying the principle of harmonic balance. Yun and Shnozuka [21] used 
nonlinear Kalman filtering algorithms for identification of MDOF nonlinear structural 
dynamical system. 

In practical engineering, many systems have multiple input variables and multiple 
output variables. Such a system can be said to be a multiple input and multiple output 
(MIMO) dynamical system. In this case, identification of MIMO dynamical system 
yields the learning problem of mapping between the multiple dimensional input space 
and multiple dimensional output space. (Fig. 5.1) 



X(k) 


U(k) 



Figure 5.1: MIMO dynamical system 


For identification of nonlinear dynamical system, the Hammerstein Feedback Model 
can be easily extended from SDOF case and SISO case to the MIMO case and MDOF 
case. 


5.1 The Hammerstein Feedback Model of Nonlinear MDOF Dynamic Sys- 
tem 

The HFM of a nonlinear SDOF dynamical system in discrete time domain is given by 
Eq.(4.1.24) as follow 

x(k) = — a i r(b — 1 ) — a„x(k-n) 

-+- -> ! u ( A- - 1 ) — b^ 7 x 2 (k - 1 ) 

T 111---- 
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+ b n i\u(k — n) + b b n ') p x p (k ■ — n) 

-|- bi^ 2 r 2 (fc — 1) + b b n n q i 9 (k - n) 


(5.1.1) 


If the input and output in Eq.(5.1.1) are now defined as input vector U(k ), 


U(k) = 


vi(k) 

u 2 (k) 

( k ) 


(5.1.2) 


output vector X(k) and output velocity vector X(k) 


X(k) = 


X(k) = 


x 2 (fc) 

x r (fc) 

ii(fc) 

x 2 (k) 

X r {k) 


(5.1.3) 


(5.1.4) 


The HFM of nonlinear MDOF dynamical system can now be written in matrix form 
as: 

A'(fc) = -±AiX(k-j) + ±BiV(k-j) 




j = l 


+ ^B^-i) + ^^I 3 (fc-j) 

j=i 

i=i 

+ fr'A-’(t-i) 

j=i 


(5.1.5) 
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where is a r x r parameter matrix. B{ is a r x m parameter matrix. B 2 , B 3 , 
B J p , Tl, are r x r diagonal parameter matrices ( j = 1 ■ • • n ), and 

EU 1 EU - i) 

x 2 (k-j)= : 

_ EU EU -i) 

eu eu eu - i) 

A r 3 (fc-i)= ; 

_ EU EU EU - >K(* - iK(* - i) 

’ E- • • • EU 7i,— »,*»! (*-»•♦• **(* - J) 

A p (fc-»= i 

EU • • • EU Hi l .. 4 ,*ii(k -])■■■ ~ 3) 

EU 1 EU ~ - j) 

x\k-j)= : 

EU E' 4 in(fc-i)i,(^-j) 


C..1 ••• s, dUM* i) 

X’(k-j)= : (5.1.6) 

where j = 1 • • • n. 

For nonlinear MDOF structural dynamical system, the order n is equal to 2. The 
HFM of noplinear structural dynamical system in discrete time domain is 

X(k) = -A 1 X(k-l)-A*X(k-2) + B*V(k-2) 

+ B\X 2 {k - 2) + BlX*{k - 2) + • • • 

+ B p X p (k — 2) 4- T\X 2 (k — 2) + • • • 

+ TlX q {k-2) ( 5 . 1 . 7 ) 

To estimate parameters, the HFM in discrete time domain at i th degree of freedom can 
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Figure 5.2: The Hammerstein Feedback Model of structural nonlinear MDOF dynam- 
ical system 

be denoted in following form. 

i i{k) = — A}X(k — 1) — A 7 X(k — 2) + (B{)jU{k — 2) 

+ (£ 2 )i,A 2 (A: - 2) + ••• + {B 2 p ) it X p (k - 2) 

+ (r 2 ) <t A' 2 (A- - 2) + • • • + (r 2 ) t ,A «| k - 2) (5.1.8) 

where A], A? and ( B \ ), are the i th row of matrices A\t x r), A 7 (r x r) and B 7 {r x m ) 
separately. ( B ])«, • • •. (B 2 )„, (r|)„. • • •, (r 2 )„ are the i th diagonal elements of matrices 
Bl. •••. B 2 . F;. • • •, T 2 separately. The HFM of nonlinear MDOF dynamical system 
in discrete time domain can be illustrated by figure Fig. 5. 2. 

Consider a nonlinear two degree of freedom spring-mass structural dynamical sys- 
tem with cube nonlinear stiffness to show the application of HFM of nonlinear MDOF 
dynamical system in discrete time domain (Fig. 5. 3). We have p — 3 and q = 0 for the 
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Figure 5.3: A nonlinear spring-mass structural dynamical system 


HFM of this system and 


X(k) = 

U(k - 2) = 


*1 (*) 

x 2 {k) 

t»,(*-2) 

- 2 ) 







(5.1.9) 

(5.1.10) 

(5.1.11) 

(5.1.12) 

(5.1.13) 

(5.1.14) 
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(5.1.15) 


where 


and 


where 


Bl = 


bn 0 
0 b% 


X\k-2) = 


x\{k - 2) 
xl(k - 2) 


x\(k- 2) = EE 7 ili3 x] i (k-2)x t2 (k-2) 

* 1=1 *2 = 1 

= lh{xi{k - 2)) 2 + (tJz + 72i) x i( fc - 2 )**(* ~ 2 ) 

+722( a '2(^ — 2 )) 2 


*S(*- 2 ) = Eii^(^ 2 )#- 2 ) 

* 1=1 *2 = 1 

= 7 1 2 1 x?(fe - 2) + ( 7l 2 2 + 72i)*i( fc - 2 M* - 2) 


+ 0 2 2 (x 2 (fc-2)) 2 


X 3 (Jt - 2) = 


*?(fc-2) 

xl(k - 2) 


— 2) — EE E 7, 1 1 i2t 3 I *i(k — 2) x iAb 2) x i 3 {k 2 

»1 =1 t2 = l »3 =1 

= 7m(*l(*= - 2))’ 

+(7ll2 + 7l21 + 721l)( I l(^ — 2 )) X ^(k — 2 ) 

+ ( 7 l 22 + 7212 + 722 l ) X l (^ — 2 )( X 2 — 2 )) 2 

+7222( X 2(^ — 2 )) 3 


and 


x\{k - 2) = E E E 7*w,*m(* - - 2)* i ,(fc - 2 

*1=1 *2 = 1 *3=1 

= ll 2 n (x 1 (fc-2)) 3 


(5.1.16) 


(5.1.17) 


(5.1.18) 


(5.1.19) 


(5.1.20) 
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+(7i\ a + 7i2i + 7 i ! u)(*i(* “ 2)fx 2 (k - 2) 

+ (7l22 + 7212 + 7221 )® i (^ — 2)(x 2 (fc — 2)) 

+ 7 2 22 (z 2 (* ~ 2)) 3 (5.1.21) 

Substituting A'(fc), U(k — 2), A 1 , A 2 , B\, B\, B \ , A r2 (fc — 2), X 3 (k — 2) into Eq.(5.1.7), 
we have the HFM of system in discrete time domain. 

sci(fc) = -a} u Xi(k - 1) - a\ 2 x 2 {k - 1) - a^x^k - 2) 

—a 2 2 x 2 {k — 2) + b 2 \ui(k — 2) 4- fc 12 t/ 2 (fc — 2) 

+^7i> 1 (* — 2)) 2 4- fc 22 ( 7 J 2 + 72i) x i(^ ~ 2 )i 2 (A- — 2) 
+&n722M* - 2)) 2 + b\l 7 m(*i(* - 2)) 3 
+^ii(7m + 7m + 7 2 h)( x i(^ - 2)) 2 * 2 (* - 2) 

+^n(7i22 + 7212 + 722i) 1 *i(^ — 2)(a ' 2 (k — 2)) 2 

+«7ii(**(* - 2)) 3 (5.1.22) 

x 2 (k) = — a 21 «i(fc - 1) — a\ 2 x 2 (k — 1) — a 21 Xi(fc — 2) 

—a\ 2 x 2 {k - 2) + 6 2 lu!(fc - 2) + &”u 2 (fc - 2) 

+fc 22 7 2 i( x i(fc - 2)) 2 + fc 22 ( 7l 2 2 + 72i) 2 *i(* - 2 )x 2 {k - 2) 

d"5 227 22® 2 (^ — 2) + 5 22 7m( x i(k — 2)) 

+lg(7?.. + 7?* + 7,\, )(*■(* - 2))’^(fc - 2) 

+ 1? I2 + 7fn)*i(* - *)(*>(* - 2)) ! 

+ ^22^22 2 ( I 2(^ — 2)) 3 (5.1.23) 

Let 

a l — a ll 
a 2 “ a i2 

a e = ^n(0i2 O 21 ) 
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a u — fcn(7n2 + 7i2i 72ii) 
a 12 = &ll(7l22 “I" 7212 722l) 

a 13 = ^ll7222 

and 

^1 = a 21 
62 — ^22 

b& — £> 22 (7i2 "b 721 ) 

^11 = ^ 22(7112 d" 7l21 "b 7211 ) 

fcj 2 = ^22(7i 22 “I" 7212 "b 7221 ) 

^13 = ^227222 

Then, the HFM of system can be represented by the following forms. 

x-[(k) = — aiXi(k — 1) — 02 X 2 ( 1 * — 1) — & 32 i{k — 2) 

— a 4 x 2 (k — 2) + a.sUy(k — 2) -f a 6 u 2 (k — 2) 

+a 7 (xi(k — 2)) 2 4- a»xi(k — 2)x 2 (k — 2) 4 ag(x 2 (k — 2)) 2 
+aio(ii(fc - 2)) 3 4 a n (xj{k - 2)) 7 x 2 (k - 2) 

+ai 2 xi(fc - 2)(x 2 {k - 2)) 2 + a 13 (x 2 (fc — 2)) 3 

and 

x 2 (k) = -b lX i{k - 1) - b 2 x 2 {k - l)-b a x 1 {k-2) 

—b 4 x 2 (k — 2) 4 hui(k — 2) 4 b 6 u 2 (k — 2) 

+b 7 (xi{k - 2)) 2 4 bsX-iX^k - 2 )x 2 (k - 2) 4 b 9 {x 2 {k - 2)) 
4fcio( I i(^' — 2)) 3 4 bii(x\(k — 2)) x 2 (A’ — 2) 

+bi 2 x\(k - 2 )(x 2 (fc - 2)) 2 4 b 13 (x 2 {k - 2)) 3 


(5.1.24) 


(5.1.25) 


(5.1.26) 


(5.1.27) 
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where parameters a l5 a 2 , • • •, 013 and hi, h 2 , ■ • h 13 are to be identified. 

5.2 Estimation of Parameters of HFM 

From Eq.(5.1.8), the HFM of nonlinear MDOF structural dynamical system can be 
considered to be r submodels at i th degree of freedom (t = l,2,---,r), in which the 
U(k) is the input of system, and X(k ), A (k) are the output of system. If 'hjijts’ 
... are equal to zero for *1 # i 2 , ii ? i 2 ± i 3 , * * for i th degree, the Eq.( 5.1.8) can be 
represented in following form. 

x i{k) = — a' n Zi(fc — 1) — o\ 2 x 2 (k — 1) + ' ’ * ~ a lT x r {k — 1) 

— a' 21 ii(h — 2) — • • • — a l 2T x T (k — 2) 

-f a‘ 31 Ui(fc — 2) + • • • 4- 0-3m U rn{k — 2) 

■+• a4 1 (ii(A: — 2)) 2 + ■ • • + a4 r (*r(fc — 2)) 

_L ... 

+ a\ p+2)l (^(k - 2)) p + • • • + a (p+2) M* - 2 )) P 

+ a (p+3)l(*l(k ~ 2 )) 2 + ^ a (p+3)r( I *-( fc “ 2 )) 

+ ... 

+ a’ (p+ , + i)i(ii(fc - 2)) 9 + ¥ o.\ p+q+1)r (x r {k - 2)) 9 (5.2.28) 

where i = 1, • • • ,r, p is the order of nonlinear stiffness and q is the order of nonlinear 
damping. By taking JV + n samples of the input U{k) and the outputs X{k) and X{k) 
and substituting the measured data into Eq.(5.2.28), we obtain a set of linear algebraic 
equations with unknown variables, which are system parameters. The set of equations 
can be written as following form in matrix. 

X' = A'Q' (5.2.29) 
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where i = 1, • • • ,r and 


X' = 


Xi{k) 

Xi(k + 1 ) 


Xi(k+N) 


(5.2.30) 



-x x (k - 1) 

-x 2 (fc - 1) 

-x r (fc-l) 

-*i (k) 

-x 2 (fc) 

• •• —x T {k) 

- Xl {k + N-1) -x 2 (k + N- 1) 

••• -x r (fc + AT-l) 

~x,(k- 2) 

—x 2 (k — 2) • 

— x T (k — 2) 

- 1) 

-x 2 (k - 1) ■ 

—Xr{k — 1) 

-x,(k + N-2) -x 2 (k + N- 2) • 

• ■ - Xr (k + N- 2) 

Ui(/c — 2) 

u 2 (k — 2) — 

u to (A; - 2) 

i/i(fc - 1) 

u 2 (fc - 1) • • • 

ti TO (/r - 1) 

Ul {k + N- 2) u 2 (fc+X-2) 

+ N - 2) 

(^(fc-2)) 2 

(*»<*- 2)) ! 

(Xr(k~2))‘ 


(ijft - 1)) ! 

( x T ( /r 1 ) )" 

M* + JV - 2)) 2 (x 2 (/c + N — 2)) 2 

... ( Xr (k + N-2))‘ 

(x,(*-2 ))” 

(x 2 (fc-2)) p 

(x r (k- 2)) p 

(x.(fc-i)r 

(* 2 (*-l)) p 

... (x r (fc-l)) p 

{ Xl {k + N- 2)) p (x 2 (fc + A r -2)) p 

... (*,(*• + A r -2)) ; 
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(^( fc - 2)) 2 (x 2 (k- 2)) 2 

(ii(fc-l)) 2 (x 2 (k- l)) 2 

(ii(k ■+ N — 2)) 2 (x 2 (k+N- 2» 2 

( x a ( fc - 2 )r 

- 2))« 

(ii(fe-l )) 9 

(*!(*- 1))' 

(x,(fc + AT - 2))« (i 2 (fc + N - 2))’ 


V 

•• 

1 


a !r 


a 21 

©' = 

a 2r 


a 31 

‘ 



a 41 


. a (P+9+l)r . 


(*r(k - 2)) 2 

(ir(k-l )) 2 

(ir(k + N - 2)) 2 

(ir(k- 2 ))* 

( Mfc - 1)) 9 

(i r (t + A^-2)) 9 


(5.2.31) 


(5.2.32) 


If we have only white noise, the problem of the identification of the parameter vector 
©' from Eq.(5.2.29) becomes a standard least square problem. Then, the parameter 
vector 0' can be estimated by SVD method. 


©’ = (A') + X' 


(5.2.33) 


where (i = 1,2, • • • ,r). 
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X 



Figure 5.4: A two degree of freedom spring-mass nonlinear dynamical system 
5.3 Numerical Examples 


A two degree of freedom spring-mass nonlinear dynamical system is considered in 
Fig. (5. 4). Mass and stiffness matrices are [22] 


The damping matrix is 


[M] = 


1 0 
0 1 



id 


0.1 -0.05 
-0.05 0.1 


(5.3.34) 


(5.3.35) 


(5.3.36) 


and the nonlinear vector [N and force [F are 


O.lxj -f 0.1 (j"i - r;)'' 

O.lx? + 0. l(a* 2 - xi) 3 
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cosuit 


(5.3.38) 


[*■] = 


0 


The differential equation of motion of the system is 


mi-*] + [cm + mm + w = m 


where 


Example 1: 



'i 2 



i i 
x 2 



x 2 


(5.3.39) 


(5.3.40) 


(5.3.41) 


(5.3.42) 


In the first example , damping is neglected, Eq. (5.3.39) becomes 

[M]\ X] + [A'][A r ] + [ N ) = [F] (5.3.43) 

The Eq.(5.3.43) can be rewritten as 

*i + 2xi — x 2 4- 0.2** — 0.3*i* 2 

+0.3*i*2 — 0.1*2 = cosu> t (5.3.44) 

x 2 — 2*2 + 0.2*2 — 0.3*2*i 

+0.3*22*1 - 0.1*1 = 0 (5.3.45) 


In order to create a set of simulated experimental data, Runge-Kutta method is 
used to numerically integrate the equations (5.3.44), (5.3.45) and find *i (fc), x 2 (k) with 
initial conditions *i(0) = 1, *i(0) = 0, 2*2(0) = 0, *2(0) = 0. Then, we have input 
data cosu ik and the output data *1 (k) (Fig. 5.5), x 2 {k) (Fig. 5.6) that can be used 
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for identification of the system. 


According to Eq.(5.1.26) and Eq.(5.1.27), the HFM of nonlinear dynamical system 
Eq.( 5.3.43) in discrete time domain can be assumed as following forms. 

For X \ : 

Xi(fc) = — Oi«i(fe — 1) — 02*l(fe — 2) 

+ a 2 x 2 (k — 2) + a 4 (ii(fc - 2)) 3 
+ a-s(xi(k — 2)) x 2 (k — 2) 

+ a 6 x a (fc-2)(x 2 (fc-2)) 2 

+ a 2 (x 2 (k — 2)) 3 cosu?(k — 2) (5.3.46) 

x 2 (k) — — bix 2 (k - 1) — b 2 x 2 (k - 2) 

+ b 3Xl (k - 2) + b 4 (x 2 (k - 2)) 3 
+ b 5 (x 2 (k-2)) 2 x 1 (k-2) 

+ b 6 x 2 (k - 2)(xi(k - 2)) 2 

+ b 7 {x 1 {k- 2)) 3 (5.3.47) 

where Oi, a 2 , ■ • • , a 7 , 6j, ■ • • , 67 are unknown parameters. Assume u; = 0.5, sampling 
time period At = 0.05, and 502 samples of input and output are taken. The results 
estimated by SVD method have been listed in Table 5.1, and Table 5.2. 


Example 2: 

The second example is Eq.(5.3.39). A linear damping term is considered. The 
needed experimented outputs x a (fc) (Fig. 5.7) and x 2 (k) (Fig. 5.8) are again simulated 
by using Runge-Kutta method. The HFM in discrete time domain is 

x\(k) = -ayi^k - 1 ) - a 2 Xi(k - 2 ) 
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displacement 



Figure 5.5: The response (/) of two degree of freedom mass-spring nonlinear system 
without damping 


oo 



1 

0.8 

0.6 

0 . 4 ^ 



time (sec.) 


Figure 5.6: The response x 2 (t) of two degree of freedom mass-spring nonlinear system 
without damping. 


p 


Real P 


Est. P 


Error 


«1 

-1.995 

-1.995 

0 


1 

1 

0 

*3 

1 



a 4 

-0.2 

-0.1994 

0.06% 

^5 

0.3 

0.29893 

0.35% 


-0.3 



a ~ 

0.1 

0.0997 

0.3% 


Table 5.1: Estimated parameters 


P 

Real P 

Est. P 

Error 

B 

-1.995 

-1.995 

0 

t>2 

1 

1 

0 

b 3 

1 

0.99952 

0.05% 

m 

-0.2 

-0.1992 

0.07% 

b s 

0.3 

0.298169 

0.6% 


-0.3 

-0.298756 

0.4% 

B 

0.1 

0.099566 

0.4% 


Table 5.2: Estimated parameters 
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p 

Real P 

Est. P 

Error 

a l 

-1.99 

-1.99002 

0.005% 

a 2 

0.995 

0.99501 

0.0012% 

^3 

2 

1.98684 

0.66% 

a 4 

-1 

-0.99476 

0.5% 

a S 

-0.2 

-0.198628 

0.69% 

a 6 

0.3 

0.30208 

0.69% 

a? 

-0.3 

-0.3061 

2% 

^8 

0.1 

-0.94646 

5% 


Table 5.3: Estimated Parameters 

+ a 3 x 2 (fc — 1) + a 4 x 2 (k — 2) 

+ a 5 (x a (fc - 2)) 3 + a 6 (x a (fc - l)) 2 x 2 (fc - 2) 

+ a 7 Xi(k — 2)(x 2 {k — 2)) 2 

+ a 8 (x 2 (fc — 2)) 3 -f cosa>(fc — 2) (5.3.48) 

x 2 (fc) = — biX 2 {k — 1) - b 2 x 2 (k — 2) 

■ + bsXxik - 1) -f b A Xi(k - 2) 

+ b 5 (x 2 (k — 2)) 3 + b 6 (x 2 (k - 2 )) 2 xj(fc - 2) 

+ b 7 x 2 (k - 2)(x 1 (Jfc - 2)) 2 + b 6 ( Xl (k - 2 )) 3 (5.3.49) 

where aj, • • • , a 8 , £>i, • ■ • , b 6 are unknown parameters. 402 samples of input and output 
are taken and u> = 0.5, sampling time period At = 0.05 are assumed. The results 
estimated by SVD method are shown in Tables 5.3, 5.4. 


Example 3: 
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Figure 5.8: The response x 2 of two degree of freedom nonlinear dynamical system with 



p 

Real P 

Est. P 

Error 


-1.99 

-1.99004 

0.002% 

^>2 

0.995 

0.995016 

0.0016% 


2 

1.99138 

0.43% 

b< 

-1 

-0.99566 

0.43% 

h 

-0.2 

-0.20172 

0.8% 

be 

0.3 

0.301755 

0.6% 

b t 

-0.3 

-0.30002 

0.01% 

b& 

0.1 

0.09905 

0.95% 


Table 5.4: Estimated parameters 


In this example, a Coulomb damping force is considered. This Coulomb damping 
force is assumed as a nonlinear term in Eq.(5.3.39). 


[N) = 


0.2 | ii | ii 

0.2 | x 2 | x 2 


(5.3.50) 


The equation of motion of this system is 


\M][x] + [A'][x] + [A r ] = [F] 


(5.3.51) 


The initial conditions are: 

Xi(0) = 1.0 

ii(0) = 0 
x 2 (0) = 0 

x 2 (0) = 0 (5.3.52) 

The simulated responses of displacement and velocity are obtained by using Runge- 
Kutta method from Eq.(5.3.51) and illustrated in Fig. 5.9 and Fig. 5.10. Sampling 
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p 

Real P 

Est. P 

Error 

Cu 

0.2 

0.1988 

0.6% 

C22 

0.2 

0.20004 

0.002% 


Table 5.5: Estimated parameters 

time period At = 0.05 and 502 samples of input and output are considered. The 
estimated results are listed in Table 5.5. 


Example 4: 

In practical engineering, the real damping usually is different from design damping. 
Identification of the difference is useful for analysis, design, and control. If the mass 
matrix [AT], stiffness matrix [A"], and damping matrix [C\ are known, the difference 
of damping can be estimated by using Hammerstein Feedback Model. The difference 
of damping is assumed to be [dC). We assume a linear dynamical system as following 
differential equation. 

[M)[x] + (\C) + [dC))[x] + [A'][x] = [F] (5.3.53) 


where difference of damping is assumed as 


[dC} = 


0.05 -0.005 
-0.005 0.05 


(5.3.54) 


The responses of displacement x a (t), x 2 (t) and velocity Xj(t), x 2 (f) are obtained by 
using Runge-Kutta method from Eq.(5.3.53) and shown in Fig. 5.11 and Fig. 5.12. 
The HFM of the system is assumed as 


xi (k) + aiXi(fc — 1) + a 2 Xi(k — 2) + d 3 X 2 (k — 2) + cos 0.5(/r — 2) 

= a 4 Xj(fc — 2) -t- a 5 x 2 (^ ~ 2) (5.3.55) 
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displacement and velocity 


In 

0 . 5 - 



time(sec.) 


Figure 5.9: The responses xj and Xj of two degree of freedom nonlinear dynamical 
system with Coulomb damping 
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p 

Real change 

Est. change 

Error 

den 

0.05 

0.050038 

0.008% 

dci 2 

-0.005 

-0.0050528 

1.1% 

dc 2 1 

-0.005 

-0.00496 

0.8% 

dc 22 

0.05 

0.04994 

0.12% 


Table 5.6: Estimated parameters 


and 


x 2 {k) + b 1 x 2 {k - 1) + b 2 x 2 (k - 2) + b 3 Xi(k — 1) 
= b A x 2 (k — 2) 4- b 3 ii(k — 2) 


(5.3.56) 


where aj, a 2 , a 3, fcj , b 2 , 63 are calculated from [Af], [A] and [C]. The [dC] has elements: 

a 4 


dc u — 

dc 12 = 
dc 21 = 

dc 2 2 = 


(AO 2 

(AO* 

(AO 2 

(AO 2 


(5.3.57) 


At = 0.05 and 502 samples of the input and output are considered, then the estimated 
parameters are shown in Table 5.6 . 


5.4 Sampling time period and estimate range 

In this section, the two degree of freedom nonlinear dynamical system with Coulomb 
damping, numerical example 3, is considered to examine the effect of sampling time 
period. This example has the equation of motion, Eq.(5.3.51), and initial condition 
Eq.(5.3.52). 
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displacement and velocity 



time(sec.) 


Figure 5.11: The responses x x . x, of TDOF nonlinear dynamical system 
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Figure 5.12: The responses x 2 . x : of TDOF nonlinear dynamical system 


07 


Fig. 5. 13 and Fig. 5.14 show the estimated Coulomb damping parameters of Eq.(l) 
and Eq.(2) of Eq.(5.9) vary with sampling numbers in a quarter of period, separately. 

Fig.5.15 and Fig. 5. 16 show the estimated Coulomb damping parameters of Eq.(l) 
and Eq.(2) of Eq.(5.3.51) vary with sampling numbers in one half of period, separately. 

Fig. 5. 17 and Fig. 5. 18 show the estimated Coulomb damping parameters of Eq.(l) 
and Eq.(2) of Eq. (5.3.51) vary with sampling numbers in a period, separately. 

The sampling ranges of responses, ii, *i, a? 2 , and i 2 , of system for estimation are 
denoted in Fig. 5.19 and Fig. 5. 20. 

The results denote that the nonlinear Coulomb parameters of MDOF nonlinear 
dynamical system can be estimated even if using few samples in a small estimated 
range. 


68 



parameter 
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0 . 2 - 
0 . 19 - 
0 . 18 - 
0 . 17 - 
0 . 16 - 
0 . 15 - 

0,1 4 0 20 40 60 80 100 120 140 160 

sampling numbers 



Figure 5.13: The Coulomb damping parameter of Eq.(l) varies with sampling numbers 
in a quarter of period 
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Figure 5.14: The Coulomb damping parameter of Eq.( 2) varies with sampling numbers 
in a quarter of period 
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parameter 


0.205 

0.2 

0.195 

0.19 

0.185 

0.18 

0.175 

0.17 

0.165 

0.16 



sampling numbers 


Figure 5 . 15 : The Coulomb damping parameter of Eq.(l) varies with sampling numbers 
in one half of period 
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parameter 



Figure 5.16: The Coulomb damping parameter of Eq.( 2 ) varies with sampling number 
in one half of period 




Figure 5.17: The Coulomb damping parameter of Eq.(l ) varies with sampling numbers 
in a period 




18: The Coulomb damping parameter of Eq.(2) varies with sampling numbers 




displacement and velocity 



time (sec.) 


Figure 5.19: The sampling range 
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displacement and velocity 
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